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Abstract 



by 

ROBERT L. LOLTILL, JR* 



Nujr,erical solutions are obtained to the sixth order system of 
disturbance equations describing the hydrodynamic stabilit}^' of a 
laminar, heated, flat plate water boundary layer, including all mean 
and disturbance property variations. The results are compared v?ith 
those calculated using the fourth order system of Uazzan, Okamura, 
and Smth [4] (which assumes only trieari viscosity variation). 

Both systems predict significant boundary layer stabilisation 
(increasied minimum critical R.eynolds Number, decreased disturbance 
amplification rates, etc.) with moderate heating, but display a 
maximun and subsequent decrease as the plate surface- to-^free stream 
temperature difference is further increased. Over the normal liquid 
range of water, the results of the sixth order calculations shov/ only 
a slight enhanceroent of stability over those obtained from the foiirth 
order system. Apparent gains in stability predicted v/hen only dis- 
turbance viscosity terms are included in the calculations are offset 
by the further inclusion of density fluctuation terms. 

The insensi tivity of stability characteristics to the addi- 
tional consideration of thermal disturbances is attributed to the 
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very limited region in v;hich such disturbances are important for 
fluids of large Prandtl number such as water. 
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CI{/\PTER I 



INTRODUCTION 



To gain some insight into the mechanism of laminar-turbulent 
boundary layer transition ^ considerable experimental and theoretical 
effort has been expended in studying the hydrodynamic stability of 
a laminar boundary layer. This work has revealed that alterations 
to the mean flow by pressure gradients, boundary roughness and 
compliancy, suction or blowing, heating and the like can cause 
significant changes in stability characteristics of the flov7. The 
present study focuses on the Influence of heating on the stability 
of a liquid boundary layer, and is motivated by the favorable effect 
that surface heating may have on drag reduction of submerged under- 
water vehicles. 

An indication of the qualitative effect of surface heating 
can be obtaJLned by utilizing the inviscid point-of-inf lection 
criterion (for constant density flows velocity profiles having a 
point of inflection are unstable). At the surface of a flat plate, 
the boundary layer equations for a fluid flow with a nonconstant 
viscosity reduce to 




( 1 . 1 ) 
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If the plate is heated (such that T^^>T^) , then (dX/dy)^^<0. Because 

the viscosity of a liquid decreases v/ith increasing temperature 

(dy/dT < 0), then (dy/dy) > 0, Conversely, for a cooled plate 

w 

(T < T ), (dy/dy) < 0. As both p and (du/dy) are positive for 
w ^ w w w 

either case, it follows directly that for T T , (d u/dy ) ^0, 

2 2 

Since the curvature, (d u/dy ), is vanishingly small but negative 
as y it follows that if its value at the wall is positive, there 

must be at least one point of inflexion in the boundary layer (i.c,, 
some point at which the curvature is zero) • Thus it is anticipated 
that heating should stabilize and cooling destabilize the liquid 
laminar boundary layer. This conjecture has been verified 
theoretically by a number of authors, (Note that it is just opposite 
to the conclusion reached for gases where viscosity increases with 
temperature) , 

DlPrima and Dunn [l] cite unpublished numerical results of 

McIntosh which indicate that for the U'/o-dimensional Tollmien- 

Schlichting type of instability in a heated water boundary layer, 

the Tuiniraum critical Reynolds number, Re , . , is increased 

^ min.crit . 

tenfold over the isothermal case for a 50°C temperature difference 
(T - 10°C), 

Hauptmann [2] utilizes integral mell’ods and a perturbation 
procedure to correlate variation in a velocity profile shape factor, 
arising from s mall changes in a temperature dependent viscosity, 
with the stability characteristics for a class of single parameter 
velocity profiles. His results also predict very strong stabilization 



in water for relatively small heatings and are quite compatible 
with those of McIntosh cited above. 

The most extensive study of the stability of a heated and 
cooled v?ater boundary layer with a Falkner-Skan flov; has been con- 
ducted by Wazsan, Okainura, and Smith [3,4,5] and Wazaan, Keltner, 
Okamura, and Smith [C]* For the disturbance equation, assuming 
explicit fluid property variation in the mean viscosity and 

no fluctuating temperature field, they obtain for all their 
analyses a modified Orr-Sommerf eld equation 

(u-cX ~ - u"(X = 

+ Jj" 0 ' >- dj)j 

where ’ here denotes differentiation with respect to y/o* In con- 
trast, the variable mean-flow coefficients are obtained by sdlving 
the coupled mean momentuin and energy boundary layer equations in 
which all fluid properties are assumed to vary [?]• In general, 
their results show that while cooling the wall does destabilize 
the flow, moderate heating very strongly stabilizes it.^ Hov;ever, 
as heating is increased, R.e , . reaches a maximum with T and 

then decreases c They explain this behavior as follows [s]: 



VThilc the results of reference [s] arc qualitatively correct, 
U’ney arc quantitatively inaccurate due to an error in tlie numerical 
integration of equation (1*2) and si;hscqiicntly have been corrected 
[a, 5]. The first reference still provides the best explanation of 
the numerical procechirc used. 



"As T is incxreased above X , the velocity profile becomes 

I^T CO ' * 

more stable., whereas the variable viscosity term, y ( n), in the 
modified 0-S equation tends to destabili?.e the flow. However, at 
moderated rates of heating, the effects of the mean velocity pro- 
files u (n) and u" (n) upon Re . . are the dominant ones. As 

min . crit . 

the rate of heating increases, the negative effect of y (n) upon 

Re . increases whereas the effect of heating upon u (n) and 

mln.crit. ^ ^ 

u" (ri) and, in turn, their respective stabilizing effects on 

Re . . begin to level off. Re . _ then reaches a maximum 

mxn.crit. ^ min.crit. 

at (T ) . where the effects of y (q) and u (q) and u" (q) upon 

w crit. n ^ / r 

Re . . balance out. As T is increased beyond (T ) _ , the 

ruin. crit. w V7 crxt. 

destabilizing effect of y (q), becomes dominant and c'^it 

decreases upon increasing beyond 

Similar results and conclusions are obtained by Potter and 

Graber [s] for the stability of Plane Poissuille v?ater flov7 with 

heat transfer. Using for their analysis the same equation as 

Vazzan, et.al. and neglecting all property variation other than 

viscosity in the solution of the mean equations, they found that 

omission of the mean viscosity gradient terms from equation (1.2) 

indicated increased stability with heating. In contrast, inclusion 

of these terms show destabilization with increased heating (i.e., 

a decrease of Rc ^ )• 

min.crit . 

Tlic motivation for the present investigation stems from an 
attempt to more fully explain the stability anomalies found by 
Wazzan, Okanura, and Smith through a more complete reformulation 
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of the problem. Specifically, it is unclear why those authors 
solve the coupled mean flow system of momentum and energy equations 
to provide the variable coefficients for the disturbance equation 
(1,2), and yet neglect a disturbance energy equation coupled with 
the momentum equation through fluctuating fluid property variation. 
As an a priori estimate of the relative magnitude of all disturbance 
quantities and their derivatives is impossible throughout the 
entire boundary layer, the omission of temperature fluctuations in 
previously cited investigations is seemingly arbitrary and v/ithout 
justification in a situation involving heat transfer. Accordingly, 
the complete linearized, parallel flow disturbance equations (ex- 
cluding only negligible expansion v/ork and viscous dissipation) 
eventually will be formulated and solved, 

A simple Illustration substantiates the proceeding conten- 
tion. Assuming for the moment that the only property variation of 
importance is that of viscosity, it v^ill be shown that the stability 
equations can be written more completely as 






( 1 . 3 ) 



6 



coupled v/ith a disturbance energy equation 



t(u-c) + H ^ 



_J 

c\' Re R- 




- o< T 



(1.4) 



through the relation between the viscosity fluctuations, m, and 
the terapsrature fluctuations r. Assuming further that disturbance 
quantities, m and as well as their derivatives, m" and are 
each of comparable magnitude at some point in the boundary layer 
and realizing that u is of the same magnitude as then terms 
arising due to the presence of the disturbance viscosity such as 
m”u and c< tnu are likewise of comparable magnitude to terms 
retained by Wazzan, et.al.such as and ^ and so should 

not he discounted in the analysis* Since VJazzan, et.al. attribute 
much of the stabilization to the mean viscosity gradient terms [s], 
it would seem as if omission of these terms is ’^not v;ell justified” 
and it is not permissible categorically ”to assume that T = 0 
throughout the boundary layer” [3], 

Thus the focus of the investigation is upon establishing the 
difference between the fourth order system of Wazzan, et.al*, 
equation (1.2), and a more complete version of the sixth order 
system represented by equations (1.3) and (1.4) in which all 
property variation, mean and disturbance (acceptable within the 
locally parallel flow and boundary layer assumptions) , is 



considered . 



A complete discussiort of the assumptions made in formulating 
the problem considered in this investigation is given in Chapter 
II. The assumed fluid propert^y’-temperature relationships, numerical 
techniques, and assessment of the numerical errors incurred in the 
process are discussed in Chapter III. The results of this analysis 
are presented in Chapter IV, and finally the conclusions dra^m from 
these results are discussed in Chapter V. The Appendices are in- 
cluded to present relevant but nonessential material that would 
otherwise obscure the primary thrust of this investigation in the 



text. 



CHAI^TER II 



FORMULATION OF THE STABILITY' PROBLEM 

This investigation examines the spatial stability to small- 

scale disturbances of the steady, laminar, boundary layer flow of 

a viscous, heat-conducting liquid. Since formulation of tlie 

problem can easily be made with far greater generality than is 

eventually required, such a procedure is followed to demonsti^ate 

that V7ithiu the assumptions to be made, the additional complexities 

fail to significantly filter the resulting stability equations. 

Thus while the problem is posed for general Falkner-Skan flows 

(u = V7ith constant cross flows (v; = constant) and arbitrary 

e e 

disturbance orientation, solutions discussed in Chapters VI and V 
are specifically for two-dimensional, flat plate boundary layers 
with two-dimensional disturbances. The Falkner-Skaii flovjs are 
often referred to as ’’wedge flows” since for M > 0 they correspond 
to a uniform flov; over a yav;ed, rigid, infinite w^edge \7ith included 
angle [9, p 14l] at zero angle of attack in the neighbor- 

hood of the stagnation point. 

2 . 1 Basic Afjs umpt ions and Gov e rning Equations 

To describe this flov? field mathematically, certain basic 
assumptions are made concerning the nature of tlui fl.uid Involved. 
( 1 ) The characteristic scales for mean motion and 



,8 
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disturbances to it must be significantly greater than those 

associated v;ith molecular motions in the fluid (continuum 

hypothesis). This permits ascription of meaningful values to 

■fluid and flow parameters ”at a point’*. Further, for disturbance 

periods (2n/a0 much greater than the molecular time scale, say 

the fluid’s structural relaxation time (T ) (i.e. wT << 1), the 

s s ’ 

flow can be considered to be in statistical equilibrium so that 
molecular transport effects can be represented by transport co- 
efficients (vi, k, etc.). 

(2) Once the continuum hypothesis is established, to define 
the stress relations betrceen fluid elements, it is assumed that 

the fluid is Newtonian (linear stress variation with rate of strain) 
and isotropic (no intrinsically preferred elemental spatial 
orientation) • 

These two groups of assumptions form the foundation for the 
Navier-Stokes equations (2. 1-2. 3). Since there is no experimental 
evidence to suggest their inadequacy under ordinary conditions , 

they will be used without reservation in this analysis, 

(3) Assuming the fluid to be a liquid, the form for the 
state equations is defined. The relative pressure independence of 
all fluid properties permits their expression explicitly as 
functions of temperature alone. 
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Thus, ( d inen s i on a 1 i y ) 



Continuity: 



Dp 

if" = o 



(2.1) 



Momentutn: 



Energy : 



State: 



wnere 



n IfU 

V Dfc 



p6. ^ 



V D 

P> >y^' ) ^ ~ "Punct ion (-r) 

^ = material derivative = ^ nT 

Dt d>t J 

or. = stress tensor = - p^. + 2^4/ £-j -f 

~ rate-of-s train tensor = *2 { "*” ' ^ ' ^1 J 






(2.2) 



(2.3) 



(2.4) 






[3 - body force 



A - volumetric expansion coefficient = ^ 

V ^ \ ST*. 

- *'iriGchanical” thertnodynamic pressure 
difference. “ ~ ( > 4 - | t) 



f 



This last equation relating the mechanical pressure (negative 
of the mean normal stress, ^ c,„.) and thcnnodynamlc pressure is 

J ioTv 

valid because only the near past is considered (or (oT << 1 as 

s 

previously discussed [lO] )♦ Note that equations (2,2) reduce to 

X 2 

the Navier-Stokrs equations if - — == - (Stokes^ liypothesis. 
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valid only for an ideal monar.omic gas) or = 0, 

For liquid in a heated boundary layer, there is no justifi- 
cation for assuming either statement is valid in general. 

Lieberrnann [11] lias shown that for a number of liquids “ is 

y 

positive, considerably greater than unity in some cases, and widely 
different for different liquids. Furthermore, for low frequencies 
(cdT^ << 1), Narasirnliam [l2] notes that the viscosity ratio is 
independent of disturbance frequency, and almost independent of 
temperature [l3]. (For water, Pinkerton [l4] has shovm a •— 
variation of about 16% over the range 0-60^C as compared with a 
y variation of 118% over the same range - see Figures B.2 and B.6). 
Alternatively, for flows with heat transfer (particularly cases 
vxth large thermal gradients), when the density is identified as a 
function of temperature, it will not be constant in that flow 
field. In such a situation, then will vanish only at a rigid 
surface for steady flows. Thus if relative] y large amplitude 
density fluctuations were expected in a liquid for which was 
large (e.g,, for benzene, this ratio is O(IO^) [ll], then there 
could be a large disparity between mechanical and thermodynamic 
pressures. It will be sho^vm in Appendix A and Section 2.3.1 that 
for values of — ~ of order unity, the stability problem can be 
solved independently of this quantity. 

When the flov; is characterized by a single time scale, namely 
that associated with the frequency of the small, rapid perturbations 
on the basic flow, each of the flow and transport parameters can 
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be constructed as the sum of a time-'averaged mean and a fluctuating 
component , 



Equations describing the behavior of these components are developed 
in sections 2,2 and 2,3. respectively. 

^ ^ The Me a n Flow 

The mean equations for boundary layer flow over a v^edge can 
be deduced from equations (2. 1-2.3) (see Appendix A) after choosing 
the orthogonal coordinate system shown in figure 2.1 and making the 
following simplifying assumptions: 

(1) The basic flow is 53teady but not necessarily spatially 
homogeneous (cqn. 2.6), so the mean equations evolve from a time- 
average of the instantaneous equations. 

(2) In the momentuni equations for motion parallel to the 
surface, viscous and inertial forces are of the same order of 
magnitude. Using the scaling resulting from this requirement and 
eliminating tcxns of 0(Re^ ^) or smaller, boundary layer equations 
are obtained. 

(3) The wedge is “infinite” in extent in the z direction, so 
all mean flov7 variation is independent of the z coordinate. 

(4) Gravity is the only body force present, acting in a 




( 2 . 6 ) 



plane normal to the z axis at an angle x to the y axis. 



Figure 2.1 Coordinate system for the stability analysis of a general three-dimensional boundary 
layer flow. 
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Additional simplifications leading to exclusion of specific 
effects on the flov? (e«g*, dissipation, compression energy, buoyancy) 
will be made after examining relevant nondimensional coefficients. 

2.2.1 Mean Flow Equations 

Using the assumptions discussed in section 2.2, equations 
(2. 1-2. 3) reduce to the (time-averaged) mean boundary layer 
equations (dimensionally) 
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Before solving, several modif ications are made to the x 



momentum equation. Evaluation of this equation at the boundary 
layer edge reveals the nature of the mean pressure variation im- 
pressed on the boundary layer (dimensionally) 



-f 



5in ^ 



( 2 . 12 ) 




^3 



Further, due to the state equation (2.4), the density can be ex- 
panded in a Taylor series in temperature, so dimensionally 









or nondimensionally ^ using the volumetric expansion coefficient, 
6, defined in equations (2.5), 




T 

For small values of 3 T and «> - 1 (such as is the case for 

T 

CO 

V7ater with moderate temperature differences) , this can x<rcll be 
approximated by 



I - ^(T) = (T- 1) «.13) 



Thus, with the modifications indicated by equations (2.12) 
and (2.13), the x momentv>.m equation finally becomes (dimensionally) 




-h 



(2.8a) 
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Since all fluid properties are assuiried to be variable, 
particularly density, a modified Howarth*-Dorotnitsyn transformation 
is suggested [7], reducing the boundary layer equations to a 
nearly incompressible form. The required independent variables 
are 



Satisfaction of the continuity equation is accomplished by intro- 
ducing a dimensional stream function. such that 



X - X 






and a nondin;ensional stream function F so that 





( 2 . 15 ) 



wliere the two are related by 




Defining the displacement thickness 




( 2 . 16 ) 
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or in transformed variables 




(2.16a) 



then 




(2.17) 



T7here 






(2.18) 



a constant for specified wall and free stream temperature and wedge 



angles. 



Additionally, defining 
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whera 



Re = 
^ = 
E = 
Pr = 
M = 



Reynolds Number = 
Grashof Number = 
Eckert Number = 



6/c ^5 



Cp T 

rco CO 

Prandtl Number = 

X dUc 

Ue dx 









= constant since u is of the 

e 

M 

fom u <=: X for the wedge flov7. 
e 



At this point an assessment is made -of the various force and 
energy terms of the above equations to determine their influence 
in establishing the mean velocity and temperature boundary^ layer 
profiles, V7ith the following conclusions. 

Gr 

(1) Buoyancy effects can be neglected if — — sin x 
This is easily accomplished ana3.y tically by requiring 

(a) the surface normal be parallel with the gravity 
vector so that x “ 0> or 

(b) for large x angles (i.e. sin x ~ 0(1)), be 

Gr 

sufficiently large that — ■ << 1, since this 

Ke 

2 

parametric group varies inversely with u^. Thus, 
except for surface orientations nearly parallel to 
the gravity vector and low speeds, free convection 
can be neglected. 

(2) Dissipation and expansion (or pressure) work are 
negligible, requiring that Pr E << 1. (An estimate of these terms 
for a water flow indicate that they can never be significant while 
the boundary layer approximation remains valid). 
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As a consequence of these assumptions, the boundary layer 
equations (2.19) are explicitly independent of x. The original 
set of partial differential equations has been transformed to one 
of ordinary differential equations, so the velocity and temperature 
profiles are given by similarity solutions. 



Note that the z luoraentuin equation is not coupled v;ith the 
other two; that is, F and H are determined independently of G 
provided p, p, and k are also independent of G, 

2.2.2 Mean Flow Boundary Conditions 

The boundary layer equations in section 2,2.1 are 
supplemented with the follov/ing boundary conditions: 



The mean flov; equations can be written finally as 





( 2 . 20 ) 





(no slip at the v/all) 
(impermeab]e v/all surface) 
(constant vjall temperature) 



( 2 . 21 ) 
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as y ^ “ 



u (x) 
c 



W 

e 

xVr 



Note that although may vary with x due to a non;^ero 
v;edge angle, the free stream temperature (identical with that at 
the boundary layer edge) is assumed constant. 

In terms of the transformed variables, these conditions 

become 



at n 0 



as ri CO 



F' = G = 0 
F - 0 
H = 1 

F’ 1 
G -V 1 
H 0 



(2.21a) 



2 . 3 The Disturbance Flow 

The disturbance equations are obtained as the difference be- 
tv:een the instantaneous and mean equations (see Appendix A) . 
Simplification is accomplished by the following considerations. 

(1) Consistent with the assumptions of section 2.2, 
"infinitesimal" disturbance amplitudes permit linearization of the 
disturbance equations. (Mote that all terms of this order vanished 
in taking the mean). This of course eliminates terras with 
correlated disturbance quantities (as for the moan equations) and 
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reduces the disturbance equations to a set of linear, partial 
differential equations vrith nonconstant mean flow coefficients. 

(2) The mean boundary layer flow is assumed to be locally 
parallel, so that for a specified distance from the leading edge, 

X, these coefficients of the disturbance quantities are functions 
only of the distance normal to the solid surface, y. That is, 
the stability of a velocity profile at a given x is to be con- 
sidered independently of the rest of the boundary layer. The 
requirements necessary to make such an assumption are elucidated 
in Appendix A and briefly sumtnarised below, the amplitude of 

all disturbances and derivatives of disturbances are equal (though 
not necessarily equal to each other) , comparison of the mean flow 
coefficients in a boundary layer sense indicates that v << u , 

w and . Furthermore, if disturbance v;ave lengths 

3x 3y 

are small compared v;ith lengths characterizing changes with x of 
mean quantities, then these quantities can be considered to be 
slowly varying functions of x compared with the fluctuations. 

Again, the disturbance equations are formulated as generally 
as possible vrithin the framev;ork of assumptions specified above. 
Elimination of other terms by examining the magnitudes of their 
dimensionless coefficients vjill lead to additional simplifications. 

2.3.1 D isturbance Equations 

With the assumptions discussed in section 2.3 the disturbance 
equations for this problem become (dimensionally) 
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(2.23) 



(2.24) 



(2.25) 
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Energy: 




Notice that by virtue of the parallel flow assumptions, the 
mean flow coefficients of all disturbance terms are functions of 
y only, so these equations are now separable in terms of solutions 
of the form 

A i CoS 6 -h £^ 0/0 0 - i'") 

Q -- pe 



Noting that the assuraed disturbance form is complex (with complex 

amplitude function q) \-;hile the disturbance equations (2, 22-2, 26) 

are not, physical quantities can later be recovered by taking the 

real part of the complex solution. Thus, the disturbances are 

assumed to be sinusoidal waves of frequency w = c propagating 

at an angle 0 to the x axis. For this spatial analysis, m is 

* 

assumed real and complex to admit solutions which are temporally 

>v *;c 

only oscillatory but spatially amplified (c<^ < 0) , neutral (o(^ = 0) , 
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or damped (o<^ > 0) « 

Such a formulation v/ill reduce the set of disturbance 
equations from partial differential equations to ordinary ones in 



Thus introducing 
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(2.27) 



Sio 0 “ C 1 



v;here 



^ - Ug (x) CoS Q -h We 5(C) 0 



(2.28) 



and defining for the mean and disturbance velocities parallel to 
the surface, respectively, 



-Tcy) -f /i (y) iao<9 

/ + iQnl]Uan6 



(2.29) 



and 



V7 = 



C/ (y) t \M (y) tfln'I.p 6ao 0 
/ -/• tan tan Q 



(2.30) 
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V 



where tan Y = 



and 6 is the angle between the x axis and the 



direction of v;ave propagation, then the disturbance equations 
become: 



Continuity 



cr ( \A/-c) 



+ L 



J' + (^0) = 0 



( 2 . 31 ) 



Momentum 
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Energy 




(2.34) 



v;here differentiation is again v:ith respect to n to permit in- 
formation obtained from the mean equations to be used directly in 
the numerical solution of the disturbance equations, and v;here 
all fluid property variation is of the form 



S ~ c/ 0 4^ 
0 
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(2.35) 



The dimensionless wave number and Reynolds number are defined by 



o( “ 






and 




Vi " 



Note that the composite x-z momentum equation is obtained by 
mtiltiplying the transformed x and z momentura equations by cosO and 
sin0, and then adding. Furthermore, diJatory expansion energy 
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drops out of the energy equation with linearization ^ and so must 
be considered a second order effect for this problera. 

To simplify these equations further, the assumptions raade in 
section 2, 2.1 regarding buoyancy, dissipation and expansion energy 
are used again, with some interesting impj.ications . 

(1) Neglect of the dissipation energy terms (for \ ) 

permits the energy equation to be written as a function of and 
W rather than fj h, u, and w, thereby reducing the order of the 
system of equations by tv:o from eight to six. In other words, when 
viscous dissipation is negligible in the disturbance energy equation 
then *‘the stability of a three-dimensional boundary layer to a 
plane-wave disturbance of arbitrary orientation reduces to a 
two-dimensional stability problem governed by the boundary-layer 
profile in the direction of wave propagation and by the mean tem- 
perature profile"[15] . Thus, alterations of the external flow due 
to changes in pressure gradient or cross flow v7ould be manifest 

as changes to the variable coefficients in the disturbance equations 

(2) VThen the expansion energy is negligible (3^T^E <<1) , 
variation of the tliermodynamic pressure appears only in the momentum 
equations. Under these circumstances, explicit inclusion of the 
second viscosity coefficient, X , in the momentum equations can be 
avoided by using the ’‘mechanical** rather than thermodynamic pressure 
Such a change would not affect the mean pressure, for using 
equation (2.!i) and the ”parallcl-f lov;“ and linearity assumptions, 

p^ - p ; however, it would result in a difference between the two 
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fluctuation pressures 




or, in terms of equations (2,27) 




(2,36) 



As noted before, for the large Reynolds and small V7ave numbers 
anticipated, this pressure difference could only be significant 



On the other hand, if the expansion energy is not negligible 
(say, (P^T^)h'^O(l)) , the appearance of X can still be avoided in 
the momentum equations. Elimination of the thermodynamic pressure 
fluctuations bet\>^een the tv;o momentum equations ( 2.33) 
yields a vorticity equation ~ an Orr^-Soinmerfeld type equation 
modified by the inclusion of fluid property terms - in which terms 
containing X identically cancel [16], However, this procedure 
would require evaluation of higher order temperature derivatives 
of the fluid properties. Since the equations will ultimately be 
reduced to a system of first order equations for solution, such a 
step should be avoided if possible. 

In short, the second viscosity coefficient should influence 
only the disturbance thermodynamic pressure distribution. For 
water flow, it has already been noted that E raust be very much less 
than unity, so the first procedure for eliminating X is utilized 
for this analysis. 



for large Vcslues of xVp'"^. 
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Since the assumption was made in solving for the mean equations 
that buo 3 ‘ancy terms of order Gr/Re << 1, it is certainly reasonable 
to assume that the body force influence with terras now of order 



included, their contribution would appear indirectly through modifl-' 
cations to the mean velocity and temperature profiles. Although 
it is knovm that for flows driven primarily by body forces, heating 
from above stabili^-es v.diile heating from below destabili:^,es the 
flow through the influence of thermal instabilities, this effect 
will not be considered herein. 

Thus with the assumptions discussed above, the disturbance 
equations finally become: 
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Honentum: 
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Energy; 

Cf,[ .(w-c)-t- +e«s)H'] = 
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(2.40) 



Since these equations are not in a form suitable for 
numerical integration, they are rearranged and recast as a system 
of six, first-order equations. Denoting the dependent variables 
as 
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the six equations can be vjritten as 






(2,42) 






V7here C,, is a 6x6 coefficient matrix with the following nonzero 
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terms: 
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Z’j. is obtained from the composite x-z raomentum ecuation; 
y from the continuity equation; , from the energy equation; 

and j from the y raomentum equation, 

2.3.2 Disturbance Boundary Conditions 

The disturbance velocity boundary conditions used in this 
analysis are easily obtained by requiring that there be no relative 
motion between the fluid and the solid surface at this interface 
(no slip), and that the velocity normal to the surface must also 
vanish there (impermeability). In notational form, at y == 0, 

[j = 0 slip) 

V = u (iinpermeabili.ty) 

or at n = 0. 0^0 (2.43) 



The thermal boundary condition at the wall requires more 
careful consideration (see also discussion by Dunn and Lin [17]). 
The scaling procedure used in Appendix A indicates that for a 
region very close to the surface (y/6 << 1), the energy equation 
reduces to the one-dimensional, unsteady conduction equation 
(dimensionally) 











7y ) 



The important point to note is that longitudinal temperature 
fluctuations are negligible in this region. Thus, solving this 
equation in the solid and requiring that disturbance temperatures 
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and heat fluxes be continuous at the liquid ~ solid irite7:facc 



(y 0) > the thermal boundary condition becomes 
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and the subscript "s" 
denotes properties of the 
solid. 



It is simple enough analytically to require that the thermal 
inertial of the solid be sufficiently large that thermal fluctua- 
tions die out in the solid very close to the surface. Under these 
conditions, the appropriate boundary condition v/ould be x(o) = 0 



(i.e. 






k: _ cfkY 



^1 ) 



Realistically, however, the full 



boundary condition should be closel}’' scrutinized when applying 
this analysis to a particular solid, at a specified frequency, etc. 
For this analysis, the condition 



T(o^ = O 



(2.44) 



is chosen. 

Further it is assumed that all disturbances are bounded at 



large distances from the surface 

y — ^ CO t/ ) ^ ) 
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^ ^ Disturbance Field En ergy B alance 

Although mathematically the stability characteristics of 
the disturbance flow field can be completely described (x/ithin 
the framexv’orlc of specified assumptions) by the equations formulated 
in the three preceding sections, it is instructive also to develop 
the theory on a more physical foundation as x>?ell. Insight into the 
physical mechanisTiis at x-jork in the boundary layer is thus sought by 
examining the production, dissipation and transfer of the dis*- 
turbance energies* 

The appropriate equations are obtained by first multiplying 
each of the x, y, and z disturbance momentum equations (2.23-2,25) 
by its corresponding disturbance velocity, u, v, and X7, and then 
summing the three. After some algebraic manipulation, a meaningful 
equation is obtained and e>q3ressed (dimensionally) as: 
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v/here V ~C> 




Assumi.ng only two-dimensional mean and disturbance flovjs, 
equation (2.46) reduces to: 
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The tine-averaged version of equation (2,47) written nondimension- 
ally in terms of disturbance amplitude functions (see equations 
(2.27) is: 
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where 

r 

and * denotes the complex conjugate. 









It is instructive to examine each term at this point to 
identify the physical process it represents and assess its relative 
importance. Thus, the indicated terms represent: 

(1) the rate of change in kinetic energy of the disturbance 
field. Note that for a neutral disturbance, the averaged value 
is zero so there is no net energy change in the boundary layer. 
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(2) the energy drawn fror.i the mean motion by tlie working of 
Reynolds stresses against the gradients of the mean velocity. As 
this is the only source of energy for the disturbance field of 
the unheated flow, it can also be expected to he significant when 
there is heat transfer. Furthermore, it is expected to be greatest 
and positive in the region roughly between the v;all and the 
critical point, signifying energy transfer from the mean ^ the 
disturbance flow, but small and quite possibly negative from there 
to the boundary layer edge. 

(3) the energy extracted from the mean flow due to the 
working of the disturbance viscous stresses agaiiist the mean 
velocity gradients. By and large, it is expected to have the same 
sign as term (2) and so should represent primara ly a source of 
energy for the disturbance field. It should be most significant 
i.n regions where temperature fluctuations are largest since the 
viscosity and temperature fluctuations are proportional. 

(A) the rate at which the pressure fluctuations do work. 
Physically this energy is produced by two different mechanisms: 

(a) the periodic expansion and dilation of each fluid element 
(caused by density fluctuations) in the presence of a distux^bance 
pressure field, and (b) the movement of this same fluid element 
along a mean density gradient also in the presence of the dis- 
turbance pressure. Note from equation (2.48) that for a neutral 
disturbance, the contribution from the first mechanism will vanish 
at both the wall and the criticaJ layer. 
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(5) change in energy resulting fi*ora a nonzero mean viscosity 
gradient* As may be seen from equation (2.48), this term is of 



order 






times that of the major production terms (2). 



Since the general 07:der of decreases and Re increases with 
moderate heating, it is expected that this term should be of 
negligible significance. 

(6) mechanical energy lost from the disturbance field by 
viscous friction. This loss is manifest physically as a heat 
source which in turn produces minute changes in the disturbance 
temperature (and so density) level. Since the term is negative 
definite, it will alv;ays represent a dissipation and should be most 
important near the v;all. 

(7) energy dissipation from thermal fluctuations. Specifi- 
cally, this is the energy lost associated v;ith the dilation of the 
fluid elements due to their expansion and contraction with changes 
in density (caused by changes in temperature). This effect should 
be small in liquid boundary layers v;here the amplitude of dis- 
turbance density fluctuations is relative] y small, even with 
significant mean thermal gradients. 

(8) rate of transfer of pressure and dilatory energy, 
vorticity, and disturbance viscous shear stress. If the disturbance 
field is neutral (o<^ ~ 0) , no energy is transferred parallel to the 
plate. Further, if this term is integrated from the v;all to some 
large distance (p «"^) , the remaining energy contributions will 
vanish identically due to the homogeneous boundary conditions. Thus, 
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the term acts only t.o transfer energy from one position in the 
boundary layer to another, and for a neutral disturbance makes no 
contribution to the net energy. 

Two examples are now examined to more clearly identify the 
tcrris of interest in this investigation. First, equation (2.48) 
is reduced to a form compatible vrith the assumptions made by 
Wazzan, et.al; namely in the disturbance equations, "p = constant 



(but variable in the mean flow) , p = y (y) , and t = 0 (so p =y = 0) . 
Thus, for a neutral disturbance, the y-integrated , nondimensional, 
energy equation is composed of only the production and dissipation 
contributions: 



(Note that the variable p still appears in this equation since 



Secondly, assvining for this analysis all mean and disturbance 
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(2.49) 



differentiation is with 




quantities may vary> the equivalent equation becomes 




where 



(q) is here slightly simplified to 



CHAPTER III 



NUIMERICAL SOLUTIOI^ OF THE STABILITY PROBLEM 

^ ^ Thermo d ynamic Propert ie s and Transport Coefficients 

Before numerical boundary layer integration can be performed 
for either the mean or disturbance equations, values for the 
thermodynamic and transport properties of the liquid and their 
required derivatives must be specified at each step. Consistent 
with assumptions made thus far that the property variations are 
independent of pressure (also for the boundary layer approximation j, 
the mean pressure is constant in the boundary layer anyway) , these 
quantities can depend on position only through values of the mean 
temperature. For the mean flow, such temperature dependence 
couples the momentum and energy equations. More specifically, for 
the stability problem, accurate information is needed not only for 
the variable properties of the mean flow, but for the disturbance 
quantities r, m, and k as v/ell (see equations (2.20 and 2.35)). 

Since the solution process for the mean flow equations 
requires changes of the temperature at a particular physical 
position both during the iterative search for unknov/n wall boundary 
conditions (discussed in section 3.2) and later v^ith the choice of 
step size for the converged solution, discrete input to a numerical 
scheme of required fluid properties and derivatives is impossible. 
That is to say, this information is available in general only by 
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interpolating tabulated e:<periTnental data at each step. It 
follows, than, that the temperature derivatives must be obtained by 
numerical differentiation of the discrete sot of data points, and 
so must be interpolated as well. 

Water is chosen as the representative liquid for numerical 
purposes to maximize applicability of the results of this investi- 
gation (in \:hich thermal fluctuations are considered) to situations 
occurring in nature, but primarily to enable comparison V7xth the 
results of the fourth order stability analysis of Wazzan, etoal, 

[4]. In general to obtain the ’^best” information available, 
values for c^, p ,y , and k are taken from recently published 
measurements and/or evaluations of preexisting data. Temperature 
derivatives of these properties are obtained directly by 
analytically differentiating the empirical property-temperature 
relationship proposed for each. Note that this procedure is in 
contradistinction to first estimating derivatives directly from 
the data (\;ith the method of cubic splines, for example) and then 
smoothing any spurious fluctuations. 

A more extensive discussion of the sources and accuracy of 
input property data, alternate means considered for curve-fitting 
tabulated data and obtaining high quality derivatives, and graphical 
presentation of the properties and temperature derivatives for 
water are given in Appendix B. 



^ ^ Equations 



The equations describing the nean t\7o-ddmensional, boundary- 
layer flov; over a heated, flat plate are obtained directly from 
equations (2.20) as 



and. 



(pS F“)' ,.1 FF" = o 
h')' + f H' = o 



(3.1) 



with boundary' conditions specified by equations (2.21) 



*^ = o ; f '= f = O ^ 1 

and (3.2) 

IT| -V 00 : F ' ^ O 



The numerical solution to these equations introduced se^'^eral 
difficulties characteristic of a boundary value problem, in 
addition to those associated with the nonlinearity and coupling 
through fluid property variation. The method of Nachtsheim and 
Swigert [18] adopted for this analysis is found to satisfactorily 
deal v;ith all difficulties, especially: determination of where 
to stop the integration for the unconverged solution, and 
approximation and refinement of initial conditions. 

First, to obtain a unique solution of these equations. 



^5 



external conditions must be ’’satisfied^* for an undetermined, but 

finite n ^ ^ 5 * Asymptotic convergence (where, F‘ 1 and 

H 0 as ri ~ becomes F' = 1 and H = 0 at n ) is assured 

max 

simply by requiring velocity and temperature gradients vanish at 

the "truncated infinity" (F" = H’ = 0 at n= n )• Note from 

max 

equations (3,1) that these additional conditions cause all higher 
derivatives to vanish there as V7ell. Hov/ever, the value of q 

max 

(or equivalently 3 that q beyond which no appreciable change is 

observed in the solution) still remains unknov/n. If the a priori 

specification of q ^ is too small, boundary layer integration is 
max 

terminated too soon and boundary conditions cannot be satisfied; 

if q is chosen too large, computation titne becomes excessive, 
max 

Secondly, numerical integration of the mean flow equations 
requires specification of as many conditions at the wall as there 
are conditions to be satisfied at the boundary layer edge. Since 
the required additional wall conditions (F^'^ and H^) are initially 
unknovm, they must be first estimated and then iteratively adjusted 
until the asymptotic conditions can be satisfied, or 



F* (F" 
w 



H 

w 



^max^ 



1 



H(F", H* ; q ) = 0 
w v; max 

with the supplemental conditions 

F"(F", H'; q ) == 0 
w v; max 



H’ (F", ir ; q ) ^ 0 
w w max 



(3.3) 



In other words, the system of equations is solved as an initial 
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value problem with the x^all gradients as the specified initial 
values. Due to the extreme sensitivity of equations (3.1) to 
values of and even fair estimates of these parameters are 
likely to cause the numerical boundary layer integration to ’’blov; 
up" prior to reaching the prescribed value of Thus no means 

would be available for refining initial estimates of the vjall 
gradients • 

The mean equations are thus solved in the following fashion. 
With an estimate of initial conditions, F" and li’ Runge-Kutta and 
Adairs-Moulton integration techniques are applied over an abbreviated 
n range. At the end of this range, the initial conditions are 
refined by determuning the least-square solution for discrepancies 
between the computed variables and the proper asymptotic values. 

When successive iterations over this interval fail to appreciably 
change the initial conditions, an error estimate E, the sum of the 
squares of the deviation of computed quantities from their 
asymptotic values, is evaluated: 

E = (1 - F')^ + 11^ + F"^ + (3.4) 

If this error is sufficiently small, then the solution is considered 

to be found. If such is not the case, then the n range is increased 

and the above procedure repeated. The boundary layer is thus 

traversed in a stepvnse fashion until errors are within acceptable 

limits. The last n range then determines n automatically. 

max 

In summary, this method yielded relatively rapid, unicjue 
convergence to the proper solution and seemed to be insensitive to 
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initial estimates of the wall gradients. 

Once solution of the mean flow equations is accomplished, 
the. variable coefficients for the disturbance equations can be 
computed directly and the latter solved in the manner described in 
Che next section. However, other incidental calculations can 
easily be performed as veil. For example, to find n-* for the 
unheated case (or more specifically, vhen density is considered 
constant), from equation (2,18) 

f / rr”/ \ 

Is- = J r 1 - ^ - J (7 - r Vy = ^ 

^ 0 

For the heated case, though, the solution is not immediately 
available from variables computed in the boundary layer integration: 




(3.4) 



Instead, using stored values of p (p) and the composite Simpson’s 

Rule formula [19, p.79], the required integral is determined 

numerically. Compvited values of are tabulated in Tabic E.3 and 

plotted in Figure C.2, An analysis in section E.l indicates an 

p = 10 to be adequate for the boundary layer integration of the 
nax ^ y o 

wean equations. 
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Conversion back to physical coordinates must be similarly 
computed. From equations (2.16a) and (2.17) 



X. = 1 



75^ 



^7 






(3.5) 



SO that 



1 

<5^ 



1 



■1^ 






(3.6) 



'0 6 

a constant, once the mean flow equations are solved and the 

■ik 

boundary layer edge is specified. A plot of the )L:y\ and 



1 



relationships for the iTiean flows analysed is given in Figure C.l. 



3 . 3 Disturbance Equati ons 

The equations describing the hydrodynamic spatial stability 
of tvjo-dimensionai disturbances in a two-dimensional, heated, 
liquid boundary layer are obtained from equations (2 .37-2 .40) by 
setting 0-0 and if) = 0 (see. equations (2.29 and 2.30)), and 
requiring that Re and to = o<c be real and c< be complex. Thus, 



Continuity: 



Momentum: 



L _ r ( g c) ^ if ^ = O 






io/ (0 0 
^2 



Lof 



^ I e®J| 

X^<a)(u-c) = - ^rr'-t I 

2a. . 



(3.6) 



— 




Energy; 



/j9 







which satisfy a set of hoirogenous boundary conditions obtained from 
equations (2 . 43-2 , 45) , 




f ^ = r = o 



and 



(3.7) 



1 



CO : 




O 



Since these equations and boundary conditions are insufficient to 
establish a nontrivial solution ^ the system is cast as an 
eigenvalue problem (i.e., nonzero solutions satisfying boundary 
conditions are obtainable only for selected combinations of the 
parameters a, to, and Re). 

Since the complexity of the equations precludes exact 
analytical solutions vjithin the boundary layer, numerical techniques 
are applied. Of numerous possible methods previously utilized for 
solving this type of problem, the orthonormalized, differential 
method used for this investigation seems to be one of the most 
efficient and general [20]. Schematically, the employed scheme lias 
the following steps: 
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(1) estimate an eigenvalue; 

(2) starting from solutions valid in external region, 
integrate numerically to the wall and try to satisfy 
the specified boundary conditions there as VJell; 

( 3 ) refine the initial eigenvalue estimate and reinte- 
grate as in (2) . 

The final two steps are repeated until "satisfaction^’ of wall 

boundary conditions is attained. Although this procedure is 

conceptually quite straightforward, its implementation is not. The 

ensuing discussion elaborates on the difficulties associated with 

each of the steps enumerated above. 

As noted previously, there are four parameters associated 

v;ith this eigenvalue problem, , o< , o) , and Re« In subsequent 

R 1 

discussion, "the eigenvalue" should be Interpreted as those tv7o 
real parameters V7hich are iteratively adjusted to satisfy the real 
and imaginary parts of a wall boundary condition. The. choice of 
the tV7o fixed parameters is determined by inf onr.ation sought from 
the analysis; for example, to obtain curves of constant amplifi- 
cation, it is found convenient to fix and Re, and then search 

for the proper values of c< and oj . Although inclusion of a 

R 

fluctuating temperature field for this investigation is expected 
to some\;hat modify the existing results of V/azzan, Okarnura, and 
Smith [ 4 ] (in v;hich no such field was considered) , the initial 
eigenvalue estimate for a specified and is taken from the 
latter. Once several "exact" eigenvalues are determined, these 
can be extrapolated to estimate others in the same neighborhood. 
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VJith ail pertinent parameters specified, the eigenvalue 
problem is now converted to a boundary value problem. Thus it 
remains to mathematically describe the asymptotic condition 
specified for the boundary layer edge. Since the mean flow co- 
efficients in the disturbance equations are all assumed to be 
constant outside the boundary layer, an exact solution can be 
obtained for n 1 Specifically, equations (2.42) reduce to 





/ 



z 



3 



4 




Zs - [ 4c(ReO-c) + 



e 



or after some manipulation 




0 



(3.9) 



and 



z; - z, • o 



( 3 . 10 ) 
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where 

Y ^ = o/ 4 c c? Re ( 1 - c) 

(3^ = Re 1r ( - c.) (3.3.1 

From Eqn. (3*10) it is easily" shown that to satisfy the external 
boundary conditions (equation (3.7) ), for the j solution vector 



Z/J’ = 



for q > n cind real 3 > 0 

® (3.12) 



Thus equation (3.9) is an ordinary, nonhoraogeneous differential 
equation with constant coefficients which when satisfying the same 
conditions of boundedness has the solution 



yep ^ (> (p^-^(r'le) ^ (;;(]) g-U^-r)e) 



(3.13) 



for n > and real a, y > 0. 



From equation (3.8) and (3.13) 









for n > and real a, y > 0* 



(3.14) 
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The values for Z^, and Z^ can be expressed in terms of Z^, Z^, 

and Z^ as 




Re 







A - 

i-o! h 6;- , 

5 (3-4-a, 








-C5 + ir)Z/j’ + 



o/ cj/ 4^ y ) 




( 3 . 15 ) 




As indicated > the three acceptable, linearly independent 
fundamental solutions ^ ^ ^ ^ and C ^ i ‘ 

can be combined to obtain a general solution to the simplified 
disturbance equations outside the boundary layer. 

With an eigenvalue estimate and specified values for the 
thi'ee linearly independent solution vectors at the boundary 



(j) ^ 



6 . 



layer edge (obtained from eqns. ( 3 . 12 - 3 . 15 ) by setting 
i, j = 1 , 2 , 3 ), a fourth order Runge-Kutta method is implemented 
to numerically integrate simultaueousJy the disturbance equation 
solution vectors to the v;all, n = 0 . After each integration step, 
hov7ever, the original linear independence is found to be so greatly 
diminished that the solution vectors become computationally 
dependent before reaching the wall. Thus, insufficient information 
is available tliere to satisfy boundary conditions. 
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In complex vector space, this can be viev;ed as a rotation of one or 
more of the solution vectors from an orthogonal to a colinear 
configuration, or a reduction of the included angle from ninety to 
•zero degrees. The difficulty is explained as follows. Although 
the exponentially grov;ing solutions to equations (3.16) can be 
excluded from the boundary conditions in a mathematically exact 
way, their influence cannot be eliminated from the numerical 
integration inside the boundary^ layer. Here, truncation or round-' 
off errors provide suitable initial conditions for subsequent 
integration steps. As the computer vrill allow the most general 
solution, the rapidly grov;ing portion eventually dominates the 
three solution vectors (parasitic error). To avoid this difficulty, 
an orthoromal basis is reestablished whenever complete loss of 
linear independence is imminent. Gersting and Jankowski [20,21] 
contrast this "near-orthonormalized integration” to the ”fully- 
orthonormalized integration” used by Wazzan, et. al. [3], in which 
solution vectors are kept orthonormal at each step. 



Kutta scheme v;ith a fixed step size from the boundary layer edge 
toward the wall until, at mesh point the attendent parasitic 
error reduces the angle between any tv;o solution vectors below a 
specified limit, or 



Integration thus proceeds in the following manner. The 
three basic solution vectors are integrated using a Runge- 




^ (3.16) 



ii-1 
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where tepi'esencs the complex .inner product of the 

vectors i and j. There, the Gram-Schmidt algorithm for 
orthonorinalizing a set of vectors is applied [22], ^’Jhen 
eigenfunctions are to be recovered, the orthonormalising matrix, 
available as a nonessential by-product of the orthonormalization 
process, is also computed and stored for later use. The integra- 
tion is then continued from using as initial conditions the 
orthonormalized basis vectors, until at mesh point the 

angle betvzeen any tx;o solution vectors is again less than Q. This 
procedure is repeated until the V7all is reached, at v/hich point 
a final orthenermalization is performed. 

An evaluation of the proper step size, h^, and orthonormali- 
zation angle criterion, Q, required for this integration is given 
in sections E.2 and E.3 of Appendix E respectively. Based on that 
evaluation, the values of h^ < .02 and Q - 45^* were^ used to generate 
the results presented in Chapter IV* (Since the integration scheme 
used requires computation of derivatives at half-step intervals, 
the mean equations thus must be integrated v/itli step size - 
1/2 h^). 

At the wall, the three independent solution vectors are 
linearly combined in an attempt to satisfy the homogeneous boundary 
conditions, equations (3.7), for that particular eigenvalue estimate. 
In general, such satisfaction is not achieved iimmediatciy unless the 
eigenvaJue is ’’exact,** and so some iterative correction scheme must 
be implemented. Thus, for the vector components 



56 



Z,(o) =^CjZ/'’(o) 

j=' 

it is required that (see also [23]) 

Z^(o) 4- Z, (o) =0 



and 



(3.18) 



^3(0) = O 

where specification of the arbitrary constant ~ u ^(0) is found 
to give good eigenvalue convergence using the search procedures 
discussed below. Following Wazsan, et, al, [23], the quantity 
actually used to clearly distinguish an eigensolution from all 
others is the complex, surface-normal, boundaiy layer admittance 



X (co.c.Re) = 

(<>) 



(3.19) 



for, the smaller the "zero" components of the eigensolution 

(Z^(o), 2 ^( 0 ), and Z^(o)) are than the nonzero ones (Z^^(o) ,Z^(o) , and 

2,-(o)), the more accurate will be the eigenvalue [24, p. 280]. 
o 

Requiring both the real and imaginary parts of vanish cstablislies 
tv;o real equations in- the two adjustable parameters (or the complex 
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‘^eigenvalue’*) • 

For iterative refinement of the eigenvalues (i.e. , finding 
the zero(s) of > tv7o schemes were employed: the plane fit of 
WaEzan, et.al [23], and the quadratic complex curve fit of 
Muller [25], In the former, a plane of the form 

’’eigenvalue” = a + a-Y _ 4* a«Y ^ 

^ o 1 oR 2 Oi 

is fitted to three succeSvSive eigenvalue estimates and their 

corresponding computed admittances. The corrected eigenvalue 

estimate for the next iteration is then simply a • \7hen Y is 

o o 

sufficiently small, Muller’s method is used to expedite eigen- 
value convergence. This scheme fits a quadratic curve in complex 
space, also through the last three eigenvalue estimates and com- 
puted admittances, of the form 

2 

Setting Y =0, this equation can be inverted using the standard 
o 

quadratic formula to give an eigenvalue estimate for the next 
iteration. Note that both techniques have the advantage that no 
derivative of Y^ and only one evaluation of Y^ is required per 
iteration (in contrast, say, to Newton’s method), A complete 
discussion of the more subtle applications of these tv7o procedures 
for an efficient eigenvalue search is relegated to section E.4, 

Vfnen the ’'eigenvalue” is determined to v;ithin desired 
accuracy, the eigenfunctions can be recovered using the stored 
orthonornalizing matrix and solution vectors [22], For standardi- 
zation, the magnitude of these eigenfunctions arc adjusted so as to 
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make the amplitude of the pressure fluctuation equal to unity at 
the Vail. If desired, energy production, dissipation and transfer 
terms can then be evaluated directly (see section 2.4). 



CMPTER IV 



RESULTS 

The hydrodynamic stability of the boundary layer formed on 

a heated, constant-^temperature flat plate (M = 0) in two-dimensional 

water flow and with two-dimensional disturbances (ij; = 0 and 6 = 0) 

is examinad and the results presented in this section. Specific 

comparison is sought with the results of Wazzan, et.al [4] v;ho 

analyzed this same configuration, but who assumed only isothermal 

disturbances of the mean variables. 

Rather than regenerating in this investigation the stability 

characteristics for the entire range of wall temperatures considered 

by Wazzan, ct.al,, selective checks are made to compare the results 

of the two. Accordingly, since those authors found that increased 

surface heating at first stabilized and then destabilized the 

boundary layer flov/ (i.e., peaking of Re . . - that value of 

Pixn . crx L . 

Re below which disturbances are damped for all wave numbers and 
frequencies), with the extremums occurring at about - 140^F 
(T = 60^F) , the numerical calculations herein will be limited 

CO 

primarily to v/all temperatures V7ithin the normal liquid range of 
water and including that temperature range where the peaking occurs; 
specifically, T^ = 90% 150% and 200°F with T^ = 60^F. Vdien 
detailed characteristic behavior is to be illusti*ated , the case of 
T = 90®F will be used, as this is conceivably within the range of 

V7 
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temperatures for v;hich e>qperisientai comparison can be made. 

An accurate solution to the mean flow equations is essential 
for solving the disturbance equations. The procedure and results 
■for the mean flow solution are described in Section 3.2 and 
Appendix C. In the latter both detailed boundary layer distri- 
butions of all the mean properties and their derivatives required 
for the coefficients of the disturbance equations, and variations 
of surface velocity and temperature gradients with the wall 
temperature are graphically presented. 

The mean velocity and temperature boundary layer profiles to 
be examined in this section are repeated in Figure 4.1. Note that 
the temperature boundary layer is much thinner than that of the 
velocity boundary layer so thermal influences are restricted to 
regions close to the wall. This observation is related to the later 
understanding of the relative importance of thermal disturbances in 
boundary layer stability calculations. 

Keeping in mind the role that a linear stability analysis 
should play in only suggesting those conditions v/hich might enhance 
or delay the onset of laminar-turbulent transition (there is no 
definitive link between the point of instability predicted by such 
an analysis and the point of transition) , this investigation is 
restricted to regions wliich could conceivably still be laminar. 

Thus, although determination of stability characteristics for 
excessively large Re might be an interesting numerical or mathemati- 
cal exercise, the results are not of interest if transition has 
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I'j.gure A.l Moan vclocicv and tO!.perature profilos fox' vaxious T 
(T - 60 ''k). ' 

fry 
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occurred . 

Much of this investigation was conducted neglecting the 
effects of a fluctuating density field (r = 0) ; that is, only 
fluctuations of viscosity Tn(y) and thermal conductivity K(y) V7ere 
included in the disturbance Liomentur.i and energy equations, respec- 
tively. Such an assumption seemed reasonable due to the relative 
constancy of the density of water over the normal liquid range. 
Subsequent calculations with variable density nevertheless were 
made and their results are included. 

Unless otherv7ise indicated, the property-temperature varia- 
tion for V7ater is that specified in Appendix B.2. Although this 
variation differed slightly from that of Kaups and Smith [7] (used 
by Wazzan, et.al. ±n their analysis) as shown in Figure B.7, there 
is a marked effect on some of the calculated stability character- 
istics, as will be shown presently. 

^ ^ Stability Char ac teristics 

As suggested in previous discussion, for a specified set 
of mean flow parameters T^, and H) , the eigensolutions to the 

homogeneous disturbance equations and their boundar>^ conditions 
describe a surface in the four-dimensional parameter space 
(a^, , cl 3 , and Rc) . It is of interest to locate minima and 

maxima on this surface and to trace the path followed by an actual 
disturbance as it propagates through the boundary layer (e.g., to 
ascertain v;hGthGr an amplification region (a^ < 0) is traversed, 
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and if so, the rate of amplification expected) « The stability 
characteristics associated with nevitral and amplified disturbances 
are discussed below. 

4.1.1 Neutral Disturbances (a^ = 0) 

Perhaps the most obvious quantitative guide to the effect of 

heating on boundary layer stabilization is the amount v'ith which 

neutral curves are shifted with changes in the wall-to-free stream 

temperature difference. The movement of these stability loops is 

characterized by changes in the minimum critical Reynolds number. 

Re , , , the maximum wave number , and the maxiraum f re- 
run, crit.' K 

max 

quency, (above which all disturbances are damped) . As may 

max 

be seen from Figures 4.2 and 4.3, with only moderate heating, both 
the fourth order system of equations of Wazzan, et.al. [4] and the 
sixth order system of this investigation (v^ithout density fluctua- 
tions, however) predict a significant delay in the onset of 
boundary layer instability (evidenced by the increased crit ^ 

and a reduced range of frequencies and wave numbers which could 
become unstable. As the vjall-to-f ree stream temperature difference 
is further increased, however, both the fourth and sixth order 
systems consistently indicate a peaking and subsequent decrease of 

Re , . . Based on these calculations for water at the indicated 

min. crit . 

temperature levels, the following conclusions raay be dravm. 

(1) It would be impossible to achieve complete boundary 
layer stabilization by indefinitely increasing the heat rate; 



All curves are for = 60 F using the property- 
nperature relationships of Kaups, et.al.(7]. 
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Q> 




Figure 4,2 CorLpiirisor* of neutral stability curves cor^iputed with and without thermal disturbanc 
for various wall temperatures. 
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Flf^ure 4.3 Conjoarison of Rc , , computec! willi and \;ithout: 

r.ln.crit. ^ 

thermal din turbances an a function of T (fixed T == 60 F) . 

w 
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(2) v'/iiile the calculations v/hich inv^i^idc thermal disturbances 
do predict slightly larger values of Re . . , the quai.itative 

variation of this quantity v;ith T (for fixed T = 60*^F) is not 
altered. 

The values of Re . _ are easily obtained for specified 

ram. crit . . 

mean flow parameters by fixing a^. = 0 and then adjusting esti- 

i K 

mates of Re and w until boundary conditions can be satisfied. An 
example of this procedure is shown in Figure 4t4 and A. 5 for several 
different cases. 

For no other calculation is the influence of the chosen pro- 
perty-temperature relationship or assumed property fluctuations more 
pronounced tJ^an in this one. For example ^ from these tvjo figures 
it may be seen that: 

(1) use of the property variation of Kaups and Smith [7] 

effects a shift of Re . * first above (T - 90^F) 

lain, crit . v; 

and then belov; - 200"^?) that computed vrith the 
variation specified in Appendix B.2; 

(2) for either property variation, calculations Including 

thermal disturbances (but with r = 0) predict a 

Re . , roughly 10% greater than that computed 

tnin • cx’i t « 

using the fourth order equation of Wazzan, et.al. 

(equation (1.2)); 

(3) inclusion of density fluctuations in the calculations 
as well, offsets the gains in stability predicted above 
in (2) as T is increased. 

An exp.lauation for the relatively pronounced influence of 
density fluctuations on Re . . calculations is found in the 

distribution of disturbance properties in the stability equations. 
Although the density fluctuations may be comparatively small (by 
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as much as several orders of magnitude less than the viscosity 
fluctuations - see Figures 4.9 and 4,10, for example), they do 
appear in both the inviscid and viscous parts of the stability 
equations, whereas viscosity fluctuations are restricted only to 
the latter. As temperatures are increased, the disparity between 
the two disturbance quantities diminishes; eventually, disturbance 
density contributions dominate those from ail other property 
fluctuations in the calculations. 

Hovjever, as shown in Figures 4,4 and 4,5, while the fluctu- 
ating viscosity terns tend to enhance stability, the density terms 
degrade it. It is important to note that over the normal liquid 
range of water, results from the sixth order system of equations 
(even including density disturbances) predict a slightly more stable 
boundary layer than does the. fourth order system for which r (y) H 0. 

4.1.2 Amplified Disturbances (a^ < 0) 

The tendency to transition is also a function of arnplifica- 
tion rates and so attention will now be given to the amplified 
portion of the stability loop (a^ < 0) . 

A typical set of stability curves illustrating the topologi- 
cal nature of the eigenvalue variation for - 90^F, ~ 60°F, and 

M ~ 0 (r = 0) is shovm in Figure 4,6. These curves are similar in 
general shape to those for the unheated Blasius boundary layer 
although dramatically different in magnitude (see reference [26, p.90] 
for example) , 
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A much more significant: indication of boundary layer stabi- 



lization by heating than Re 



inin, crit • 



is given by the inaximuTn 



amplification rate, established for specified mean flovx 



conditions. The parametric location at which this maximum occurs 
is determined by a more involved procedure than that for finding 
Re . . , as illustrated in Figure 4,7, More specifically, the 

amplification rates for specified Re are computed and individual 
local maxima determined. Follov/xng the locus of these individual 
maxima, the largest value attaii"ied is designated the maxiiiium ampli- 
fication rate and its associated wave number (and frequency) are 
obtained graphically, RTiile the flatness of this locus curve in the 
neighborhood of allows relatively accurate evaluation 
of that quantity, it does, hov/ever, preclude similar accuracy for 

(and 03^) there. Physically this implies that there is effectively 
a band of V7ave nuTiil:)ers (and frequencies) which will receive nearly 
raaximum amplification. To complete the eigenvalue set, Re can be 
determined from Figure 4,6 by again following the locus of indivi- 
dual maxima until the required value of is reached. 

Due to the numerical effort required to compute the maximura 
amplification rate, only one set of calculations Mas made for 
T = 90®F and T = 200°F (T = 60®F) , using the presumably more 
accurate proper ty-terperature relationships of Appendix B , 2(assuminig 



and 19.5 for l;he latter can be compared with tlie correspondinr, 
values tabulated by WaJczan, et.al. [4] of 28 and 20, and the vaJue 



r = 0) . Values obtained for 




X lO' of 26.6 for the forner 



max 



maximuin amplification rate 
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for the unlieated (Blasius) boundary of 

Thus, not only is the onset of disturbance, instability delayed 
and the range of disturbance wave numbers and frequencies receiving 
amplification reduced as pointed out in the preceeding section, but 
the amplification rates to v/hich disturbances v;ill be subjected are 
also dramatically reduced. These features all point to a more 
stable boundary’' layer flow. 

4 . 2 Disturbance Amp li tude and Ener g y Distributions 

Finding no major modifications to the stability characteris- 
tics of a heated, v;ater boundary layer by including in the calcula- 
tions the effects of a disturbance temperature field, an explanation 
is sought from the boundary layer distribution of disturbance 
quantities and the energies produced or dissipated by them. 

Due to the homogeneity of the equations and boundai^^ con- 
ditions, solutions to the disturbance equations are independent of 
an absolute amplitude (ss section 2.3.1). Thus, as a comparison 
is sought bet^«;een disturbance levels for both the same and different 
mean thermal boundary conditions, it is necessary to resort to some 
physical reasoning to determine a normalizing quantity that would 
remain virtually constant in the comparison. Since there is no 
strearrwise pressure gradient (i.e., M = 0 for nil Re), tlie mean 
pressure within the boundary layer is assumed constant, and the 
disturbance pressure is nonzero and slowly varying throughout 
(see Figures 4.8 through 4.10, for example), the amplitude, of tlic 



Amplitude 



7A 





Fif.urc ^*.8 NGutrnl dinturbnnco atnplltudo profilnrs for the 
unheated (Ulasius) Itoundary layer. 



Amplitude 
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Figure A. 9 NeutraJ disturbance amplitude profiles for the heated 
boundary layer: T « 90^^F and T 60^ i'. 
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Fif.urc A. 10 llcucral dlsturbanco iir..plitudc nrofllr;; for the heated 
bouadnry layer: « 200^F and ^ 60 F, 
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disturbance pressure at the wal3 , |tt (0) j , vas chosen to normalize 
all disturbance quantities. Thus all comparisons made between 
different eigensolutions in the remainder of this section are pre- 
dicated on this assumption. 

To further investigate the importancev.-of density fluctua- 
tions, they are included in all calculations made in this section. 

Typical amplitudes of disturbance quantities found in this 
analysis are presented in Figures 4»8 through 4.10 for wall tempera- 
tures = 60°, 90°, and 200°F (T^ = 60°F) , respectively. For 
standardization, the eigenvalues associated with each set of eigen- 
functions are all located on the upper branch of their respective 
neutral stability curves at frequencies near that receiving maximum 
amplification. 

After studying these figures, the follovjing observations are 

made. 

(1) As anticipated from the relatively large Pr for water 
(i.e., such that the thermal boundary layer thickness is much less 
than the velocity one), in the temperature range considered, thermal 
effects arc restricted to a region very close to the v;all. Since 
thermal fluctuations must vanish at the wall, even v/ithin this 
region, they attain fairly large amplitudes only within a far narrov;- 
er range. 

(2) By far the most dominant property disturbance amplitude 
is that of viscosity. However, for relatively small surface heating 
(T^ = 90°F) , even the thermal conductivity disturbance amplitudes is 
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coiTjnensurate v;ich that of the normal velocity in some regions. 

For large heating (T = 200°F) , when density fluctuations are 

w 

significant, they are seen to be larger in amplitude than both the 
thermal conductivity and normal velocity disturbances* 

(3) With increased heating, the ratio of | f I / ip 

‘ ‘max ‘ ^ ^ ‘max 

increases significantly. Further, assuraing that the disturbance 

v;all pressure amplitude is relatively insensitive to changes in 

within the range 60^ ~ ~ 200®F, 1^ increases raonotonically 

V7ith T while Ip 6| exhibits a maxiinum. 

V7 I r T 

(4) Just as the stability characteristics are only slightly 
modified by the inclusion of thermal disturbances in the analysis, 
so too are the disturbance velocity and pressure amplitude pro- 
files. Even the disparity shown (except in the region where thermal 
disturbances obviously dominate) can be attributed in part to the 
difference betv/een the neutral eigenvalues required to satisfy the 
fourth and sixth order systems. 

Once the eigenvalues and eigenfunctions are obtained by the 
procedure outlines in section 2.3.1, it is a simple matter to com- 
pute the energy balance developed in section 2.4, Hov:aver, for ease 
of analysis, only those terms contributing a net increase or de- 
crease of energy of a neutral disturbance are evaluated (i.e., all 
terms of equation (2,50)). The results of sudi calculations are 
presented in Figures 4, H through 4,16. For each of the cases 
T 60®, 90®, 200®F (T - 60® F), the eigenvalues selected from the 
upper and lov/er branches of the neutral stability curves (see 
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Fif.ure A. 11 Energy d ) s crihutiojia of a noutral disturb.:' nee In thn 
unheated (Blaslus) boundary layer (lower branch). 



Amplitude 
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Fij'urc 4.12 Enerpy <ilr:t rlbuCioim uf n neutral disturbance la the 
unhoatcd (Blaslus) boundary layer (upper branch). 
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Fj^^ure A. 13 Fncrj-ry rils tribution^: of a ne^itral disturbance in the 
heated boundary layer: =* 90 F and p 00 !• (Joi;er branch). 
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Figure 4.15 Energy distri!>utions of 
heated boundary layer: T 
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'. neutral disturbance in the 
2C;0’^' and T - CO^^F (lo'/er branch). 
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T = 200 F 
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a = .20806 

K. 

u = 1.3480x10 



14383 
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Figure 4.16 Energy dist ribtitions of a neutral d is turb.'ince in tbe 
heated boundary layer: T =• 200^F and T « GO^^F (upper branch). 
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Figure 4,1) again are for frequencies for which disturbances are 
subjected to approximately Tnaximuin amplification v;hen traversing 
the loop. 

From these six figures, the following observations are made. 
The ”terms" referred to are those designated in equation (2.50). 

(1) The dominant terras in the neutral disturbance energy 

balance are the traditional ones: The Reynolds 

stress production (term (2)) in the region between 
the V7all and the neighborhood of u - Cp , balanced by 
the dissipation (term (6)) v;hich is extreraely large 
near the wall. 

(2) Dissipation due to a nonunifora mean disturbance 
density (term (7)) is always negligible. Although 
not nearly as small, energy production due to a mean 
viscosity gradient (term (5)) is also insignificant. 

Neither are large enough to appear in Figures 4.11 - 
4.16. 

(3) The disturbance shear stress (term (3)), although 

dissipative in some regions, represents primarily 

energy production for the disturbance field. In 

all the cases examined, it produced roughly 10% of 

the energy, with the Reynolds stress production 

(term (3)) supplying the rest. Recall that this 

is approximately the same amount that Re . 

^ mrn.crit. 

was increased for sixth order system calculations 

(r = 0) beyond that calculated v;ith the fourth 

order system (see Figures 4.4 and 4.5). 

(4) Pressure production (term (4)), also resulting 
from density nonunj^f orraities , primarily gives_rise 
to production for u > c^ and dissipation for u < c^^. 

Since its integrated value over the boundary layer is 
nearly zero, this term is not thought to be of 
significance in the stability of this flow. 

To check the accuracy of the computed energy production and 
dissipation terms, each was first numerically integrated from the 
wall to the boundary layer edge and then the resultant summed. With- 
in the accuracy sought, this sum was zero as required by 
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equation (2,50) • Furthermore, it was found that the percent of 
positive or negative energy contribution from the production terms 
(such as the Reynolds stress) remained virtually unchanged with 
changes in the v;all temperature level. 

It is especially important to note the extant of the region 
in which significant energy contributions take place. For the 
heated boundary layer, in contrast to the imheated one, all 
energetics are restricted to a region very close to the wall. Out- 
side of that region, energy terms tend to vanish, signifying that 
there is very little correlation betv;een disturbance quantities over 
most of the boundary layer. 



CHAPTER V 



SU>i2iARY OF R ESU LTS AND CONCLUS IONS 

This investigation has examined the sixth order system of 
equations describing the hydrodynamic stability of infinitesimal 
disturbances in a heated, flat plate, water boundary layer flow. 
Throughout y comparison has been made with the fourth order 
analysis of this same problem as previously done by Wazsan, 

Okamura, and Smith [4] . 

The time-averaged (mean) equations of laminar flow arc 
simplified by making the usual boundary layer assumptions, ai-id by 
neglecting buoyancy forces, viscous dissipation and expansion 
energy. They are then transformed from partial to ordinary 
differential equations using Ilov^arth-Dorodnitsyn type similarity 
variables and their resultant solutions used to provide the velocity, 
temperature, and other profiles required for the disturbance equa- 
tions. The latter are found as the linearized difference between 
the instantaneous and mean equation, after neglecting here as well 
buoyancy, dissipation, and expansion energy terms, assuming the flow 
to be locally parallel, and finally using a normal modes analysis 
to obtain ordinary, linear differential equations. The two sets of 
equations, mean and disturbance, thus derived differ from (hose of 
\lazzax)y et.al. by the inclusion of a disturbance energy equation 
coupied to the other disturbance ecjuations. This permits density, 
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viscosity^ and thermal conductivity fluctuation effects to be 
considered in the stability problem. 

Both the fourth order system of Wazzan, et.al. and the sixth 
order system of the present investigation predict significant 
boundary layer stabilization vdth moderate heating, but display a 
maximum and subsequent stability decrease as the plate surface-to- 
free stream temperature difference is further increased » The 
computed stability characteristics (minimum critical Reynolds Number, 
maximum amplification rate, maximum vrave number and frequency sub- 
ject to amplification, etcO obtained from each system differ only 
slightly, however. In particular, the peaking of each characteristic 
vrith V7all temperature (at about = 140^F v;ith fixed = 60®F) is 
neither eliminated nor significantly altered by the higher order 
system. Thus, complete stabilization apparently cannot be achieved 
by indefinitely increasing boundary layer heating. 

The insensitivity of these results to the additional con- 
sideration of thermal disturbances is attributed to the extremely 
limited region in which such disturbances. are important. Due to the 
relatively large Prandtl Number for water and homogeneous dis- 
turbance boundary conditions, the boundar^^ layer distributions of 
thermal disturbance amplitudes ’’spike” very close to the wall. 

In all cases considered, the sixth order stability system is 
found to predict a slightly more stable flov; than the fourth order 
system, as evidenced by larger values of the mdnii^um critical 
Reynolds Number and reduced values of the maximum amplification rate. 
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However, it is postulated that for wall temperatures beyond the 
norinal liquid range of water (with = 60®F) j, the presence of 
density fluctuations Mill eventually negate the increased stabilit}’- 
predicted with only viscosity fluctuations assumed, and actually 
degrade the degree of stability predicted by the fourth order 
system. That is, v;hen thermal disturbances are assumed, it is 
necessary to include property disturbances to obtain an accurate 
assessment of the boundary layer stability. 

The admission of thermal fluctuations also has revealed the 
existence of several energy production and dissipation terms not 
present in an energy balance equation for the fourth order system. 

Of particular note is the disturbance shear stress production, 



the energy extracted from the mean and supplied to the disturbance 
flow. In contrast, when onl}^ isothermal disturbances are assumed, 
energy is supplied completely by the Reynolds stress production. 

In summary, the fourth order system of equations used by 
VJazzan, et.al, [A] seems to adequately describe the stability of a 
heated, v;ater boundary layer flov; in the temperature ranges 
considered. This may not be true, however, for other liquids, 
particularly those with Pr < 1, v;here it is conjectured that calcula- 
tions will he much more sensitive to the influence of thermal dis- 




a source v^hich accounts for about 10% of 



turbances . 
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APPENDIX A 



Detailed D eriva tion of Me a n and Disturbance Flow Equations 

In accordance with the discussion in section 2.1, it can 
be assumed that the flo\7 field is describable by the set of 
conservation equations (2*1-2. 3). Mean and fluctuating behavior 
are deduced from these equations in the following fashion. 

The instantaneous value of each quantity can be decomposed 
into a mean and fluctuating component. 



Q. = Q< (t,) * 



(2.6) 



v?here is the time-averaged mean 








T 






c^r 



(A.l) 



and V7here the time-averaged value of any fluctuating component 
is defined to be zero. 



Q 



lihi t - 

T— CO np 



I Q(l,r)dr 



= o 



For a time average to be ruoaningful^ mean quantities must be 
time-independent as already denoted by the argument of 0^ in 
equation (A.l). The use of time- rather than space-averaged means 
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corresponds to observed boundary layer behavior v:here disturbances 
grow or decay spatially rather than temporally.^ 

ThuSj substituting equations (A-l) into equations (2.1-2. 3) 
and performing a time average, the mean equations become 
(dimensionally) 



Continuity: 



J 



(A. 2) 



Moment lira: 



^ t [ If" " f:l 

= ^ B, - ( p - 4- > J 4. 2 | yiicj 
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Energy: 
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(A. 4) 



^Eventually (see section 2.3), to reduce the disturbance 
differential ^equations from partial to ordinary ones, disturbances 
of the form e/p [f S^-ctat] will be considered. 

Clearly the time average of tills disturbance is zero, whereas a 
space average is not (except for a neutral disturbance where c<^ = 0). 
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To avoid the closure problem characteristic of turbulent 
flo;;s, it is desirable that this set of equations be solved 
independently of correlated disturbance quantities. Thus, a set of 
nev7 nondimsnsional independent variables is defined 






+ -i 






(A. 5) 



X'jhere and 1^ are length scales in the jth direction character- 
istic of mean and fluctuating flov7 behavior, respectively. For 
example, at a given position, might represent a length character- 
istic of the boundary layer development, specifically the distance 
from the body*s leading edge, v^hile 1^ would correspond to the 



disturbance wave length at that point. When characteristic length 

1 , 

scales arc comparable so that 



=; 0(1), then equation (A. 5 ) 

"j 



reduces to “ x.. Here x. must be considered a measure of both 
13 1 

large and small scale variations in the jth direction, since both 



scales are in fact Identical, or formally = X. x.. Such a 

1 3 3 

situation typically arises in this analysis where^. transverse 
boundary layer scales and 1^ are equivalent (exception to be 
noted later) and equal to the boundary layer or displacement 
thickness, for example, so that y = Y. 

If it can be argued that variations in the mean notion be 



independent of small scale disturbances in each direction (i.e. 
mean equations can be solved independently of disturbance 
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correlations) , or that the mean flow with sufficiently small 
disturbances behaves as an undisturbed flow^ then Q. must be 

X 

independent of (unless, of course, 0(1))* With this 

assumption (dimensionally) 



Qd^Xi - Q. (xp * Q,- (Xj,x,.t) 



(A. 6) 



and the mean equations become (nondimensionally ) 
Continuity: 
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Momentum: 
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^x, a>x. 
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Frora these equations, several interesting observations can 
be laade, 

(1) T\-7o reference mean quantities are left unspecified by 
any physical or geometric considerations: v^j a reference normal 
velocity; and a characteristic normal length scale * The 
former can be found from the continuity equation by requiring that 
a nontrivial tvro-dimensional solution be obtainable, so 



The latter is found by comparing the relative magnitudes of terms 
in either the momentum or energy equations. In the velocity 
boundary layer, viscous and inertial terms of the momentum equations 
ai'G assumed equal, indicating that 



Alternatively, equating conduction and convection terms of the 
energy equation in Hie thermal boundary layer 
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Since attention herein is to be focused on the stability of a flow 
field, the length scale characterizing the velocity boundary layer 
is chosen. 



(2) For the assumption that the mean behavior be independent 



must be satisfied, where the amplitude of all disturbances, have 



as small as is necessary to satisfy this criterion; experimentally , 
control of the amplitude is quite a different story. 

Thus, formally assuming that: 

(1) mean flow is steady; 

(2) viscous and inertial effects are equal in magnitude, so 
that attention is focused on the behavior of a velocity 
boundary layer, and 



of small scale variation, Q = Q(X^), to be consistent, from the 
mcmentum equations, the criterion 




been assumed equal. If 1^ represents a characteristic wave 
length, then 




vfhere 



>> 1. Analytically, the amplitude can be assumed 




y 



iOi 



(3) all disturbance amplitudes are of the same order of 
magnitude and sufficiently small that correlated 
disturbance quantities need not be considered in 
solving for mean variations, or specifically 



Qo 






0 






(4) flow quantities are independent of the z coordinate 
or the body is "infinite*' in the 2 direction, so that 



L 



X 1 



L, 



and 



B,= 0 



(5) (l?o)-^^0and = OO) but << 1, the 

last condition itap lying that the disturbance wave length 
is small compared with a length characterizing changes 
in mean quantities; 

(6) (ly/L^J =0(1) and (1^/Ly)‘^ = 

( 7 ) ^ - 0 ( 1 ) 

(8) dynamic pressure is a good measure of the system pressure, 
2 



)0 p -- p U 
O 0 0 
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the mean equations to 0(Re- ) can be written as 

t*y 

Continuity: 
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Momentum: 
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Formally, the disturbance equations are obtained as the 
difference between the instantaneous and mean equations. Per- 
forming such an operation v;ith equations (2.1-2.3) and (A.3-A.5), 
the complete set of disturbance equations can be v;rjtten non- 
dimensionally as 
Continuity: 

H- (f + S p ) A + £» p A A - ^ t ^ Cs U:) =0 

^ ^0 ^ ^ ^0 u, ^ (A. 13) 
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To make this set of equations tractable, the equations are 
linearized through the assumption of **inf initesimal*’ disturbance 



amplitudes* The linearized disturbance equations are. 
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Note that the resulting equations are those of a parallel flow. 
It can be shown by this analysis that the omitted ”uon~parallel” 



- 1/2 

terms of the boundary layer equations are of order Re. ) 



necessary to retain these terms. However, since it is expected 
that heating will increase stability [1 - 7], larger values of this 
parameter are expected than for the unheated case, and the parallel 
flow assumptions will be even more accurate for this boundary layer 
f lov; . 

From this analysis, several interesting observations may 
be made: 

(1) To the linear order considered, there is no provision 
to account for thermal fluctuations of specific heat, second 
viscosity coefficient, or coefficient of thermal expansion (however 
slight they may be) . 

(2) Additional length scales encountered in boundary layer 
stability can be easily derived from dimensional considerations. 
Consider a region very close to the boundary (viscous sublayer) with 
a characteristic velocity u^^ and length 1^^^ such that u^^ << u^ 
and 1^^ << (all other conditions remaining as previously 
specified) . If the boundary can be approached so closely that 



so that if small values of Ke^ were expected, it would be 
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>> 



, then equating viscous and inertia terms 
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of the moiuentuin equation requires that 






Incidentally, the pressure terras is also of the same order^ 
so the full laomentum equation is approximated in this region by 
Prandtl’s equation [27,p 6“9] 

^ ^ i- -L II 

<bt dx 



If another characteristic velocity, u , and length 1 are 

oc ^ yc 

u 2 ^^ 

chosen so that 1 << << — ^ — and still 1 << L , then in 

OC ty 

this region, equating the inertia and viscous terms gives 
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This region corresponds t:o the "critical layer" v/here the wave 

speed and mean velocity are equal (e.g., u = c for a two-diraensional 

K 

flow) . 

For example, specifying = 90"F, = 60°F for a fwo- 

climensional flat plate boundary layer, a solution to equations 
(A.16-A.20) with homogeneous boundary conditions is found to be 
(see Appendix E) 



Re,. = 7000 
0 

= .122347 
=0 



= 3.68283 X 10 



-6 



= .210709 



C, = 0 



where c^ = u at y = .31852 
R ^ c 



Thus for the length and velocity scales with the reference length 



L chosen as the displacement thickness, 6* 

y 



and 



-2-» 51.355>>— ^ = 4.7459 » 1 

u , u 

ob oc 




.085653 



1 

6 



.19431 
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In the viscous sublayer 0 < 
16.865 whereas frora the continuity 




equation 



u 

1 o 

so that > >. 

u 




> 599 .^ 5*7 



In the critxcal layer, y - 
u 

that 6.6978 > -2~ > 3.7104 and as 

u 

u 

critical point, ~ 4.7459. 

c 

Thus all the length and velocity scales derived from 
dimensional considerations seem to be in accord with derived 
solutions. Note that these two length scales, representing regions 
in which viscosity plays a dominant role, are also consistent with 
those derived from asymptotic analyses (see Mack [27, pp. 4-4,4*-5] 
for example) . 
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Appendix B 



Themody n andc an d Tr ansport C oe fficients and 
Their Tenp erat ure Derivative s 

^ ^ Curve-Fittinp and Derivative Deternination 

Selection of the curve-f itting technique used in this 
investigation is predicated both on the method ’ s capability of 
closely approximating tabulated data and for predicting high quality 
derivatives, with particular emphasis on the latter. In the ensuing 
discussion, no exhaustive ana3.ysis of all available techniques for 
achieving this objective is intended or attempted; remarks will be 
limited only to procedures considered. 

Of course, a given set of data can easily be **fit'^ as 
closely as desired by using a least-squares polynomial of sufficient- 
ly high degree. However, if the degree is too large, small scale 
oscillations about the primary variation V7ill be evident in the 
resulting curves due to truncation and scatter of data infonnation. 
Utilization of these curves to generate derivatives, then, will 
naturally indicate even more conspicuous but completely fictitious 
^wiggles'* in the computed derivatives. For a lovjcr degree poly- 
nomial, not only may data be poorly approximated but variations in 
higher order derivatives night not even be seen. For both the 
method of least squares and other such empirical formulations 
considered (see relationships for p and p section B,2, for example) 
it is significant to note for later reference tiiat property varia- 
tion over the entire temperature range is described by a single, 
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analytic, continuous function. Thus interpolation is not between 
”exact*‘ tabulaUGCi experimental data points, but betv;een approximate 
values (which presumably are within acceptable experimental error 
limits). Derivatives are then available immediately simply by 
analytically differentiating the indicated equation. 

As an alternative to least-«squares polynomial curve fits, the 
method of matheraatical splines can partially satisfy the previously 
established selection criteria (assuming availability of a sufficient 
quantity and quality of data information). In general, a spline 
function is composed of piecewise polynomial arcs of degree n (where 
n is a positive integer, n > 3) which are joined at prescribed points 
such that derivatives up to and including order n~l are everyv/here 
continuous. Obviously such a function provides continxAity of the 
greatest number of derivatives consistent with the use. of polynomials 
of lower degree than would be required to fit all data points 
exactly by a single polynomial. Since the curvature of a mathemati- 
cal spline is most easily controlled when n ~ 3, or a cubic spline, 
the use of simple cubics seems to be an attractive v/ay to interpolate 
(in the strictest sense) experimental data points representing some 
physical relationship. As applied to this investigation, an explicit 
formula for this “piecewise cubic” is given by 



T-1 
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where the and the J (J ^ 3) distinct temperatures over the 
normal liquid range for water 



0*C = T^< < ^3 < 



< 100 





P 



j 



and where. 



- < 



0 -for T 



(T-T^‘ U- T.T, 



The J+2 coefficients a, b, c, A-, . . ^ are then deterniined by 

X 0 *'■ i 



derivatives at the ends of the temperature interval, or [28], 
The method as suggested should v;ork perfectly if data points V7ere 
exact. Not only v7ould this data be exactly satisfied at the end- 
points of each polynomial segment, but the methods minimization of 
the curvature and relaxation of the overall requirement of analyti- 
city v/ould ensure first derivatives of high quality [29]. For 
reasons intimated above, however, it is unreasonable to expect that 
data information sliould be satisfied exactly. Consequently, the pro- 
cedure outlined above Is modified to permit departures from the 
prescribed data within specified tolerances while excluding unv/nuted 
changes in curvature, Iiuncricnlly , the smoothing is accomplivSlied by 



requiring s(T^) = j 1,. 



• > 



J and specifying first or second. 
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minimizing 



cr j-( 




v;hcre C is a positive proportionality constant and is the second 
difference at the jth point defined by 






L- -■'j- 




p - p 
"j j: 




Note that if G = 0, then only the first sum would be minimized, and 
the procedure then corresponds to the method of least squares. For 
large C, primarily the second sum would be minimized in an attempt 
to eliminate undesired fluctuations [ 28 ], Further note that in 
contrast to the first methods discussed in V7hich temperature deriva- 
tives v;ere obtained after smoothing data, this one computes accurate 
first derivatives while fitting the spline function to the specified 
data. 

For the present application, the problem is thus one of 
mathematically approximating tabulated fluid property data suf- 
ficiently well that all physicalJy realistic variation is retained, 
but spurious f lucttiations arc smoot]\cd out. Note that failure 



115 



to accomplish the latter should certainly be revealed by evaluating 
higher order derivatives. As the order of the required derivative 
is increased, the ir-oro difficult this problem is to satisfy. It is 
just this difficulty that determined the form of the disturbance 
equation adopted for this analysis (see section 2.3.1)* 

As an example of the methods discussed, Figure B.l compares 

the temperature variation of the first derivative of c obtained 

P 

from various sources by various techniques. Specifically, included 
are derivatives obtained both analytically from least-square poly- 
nomial fits of Touloukian, et.al. [30] and Kaups, et.al. [7], and 
numerically using a cubic spline fit with smoothing to the data of 
Osborne, et.al. [31]. Consistent v/ith previous discussions, the 
latter category should indicate the raore accurate general trend. 
However, it v;as not used for direct input to the numerical program; 
the following discussion explains some of the considerations dictat- 
ing this decision. 

(1) Even with smoothing, questionable oscillations in deri- 
vatives are indicated - specifically, those attributed to data un- 
certainty. In short, vrith no unusual property-temperature variation 
appearing in the normal liquid range of v;ater (e.g., numerous changes 
in slope and curvature) , data can be v^ell approximated by a least- 
square polynomial of relatively low degree. 

(2) Since the amplitude of a fluctuating fluid property at 
a given mean temperature is assumed direct] y proportional to the 
first temperature derivative there (see cqu. (2.35)), it is essential 
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Fipuro B,1 FirnC tcir.poratut o d(^ri votive oi. c 



a function of 



Lenporaturc as obtained u3ia;^ different procedures. 
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t:o know v.iiether and v;here erratic behavior in this derivative 
occurs. Only in this v;ay can logical explanations be advanced for 
aberrations in eigenfunctions, Reynolds stresses, dissipation, etc, 
(3) For simplicity of evaluation and conservation of data 
storage, it is more efficient to use a lov: degree least-square 
polynomial or empirical curve fit. 

Similar analyses v/ere conducted for the three other fluid 
properties, p, k, and p. The results suggested the property-tem- 
perature relationships for use in this stability study and are 
discussed in section B.2, 

^ ^ Property-Temperature V ariation for V/ater 

Since a major contention of this study is that a fluid's 
property-temperature variation plays an important part in establish- 
ing its flow stability characteristics in a heated boundary layer, 
accurate values for thermodynamic and transport properties is con- 
sidered essential. Tlie literature survey prompted by a search for 
the ’'best” available property information for water led to the 
following results and conclusions. 

(1) Rather than digressing on an extensive evaluation of 
data (v;hich would entail an analysis of experimental tedmiques 
and the quality of information available from each, assessment of 
systematic error limits in reported data, correlation of the infor- 
mation judged most reliable, etc.), recourse is made to recent 
studies of such a nature. Although, as evidencecl by Figure B.2, by 
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far the largest property variation with temperature is by the vis- 
cosity, equal care is exercised in specifying the behavior* of each 
of the properties. In this vray, discrepancies betv;een the results 
obtained in this numerical study and those occurring in nature are 
less likely to be attributed to inaccurate property information. 

It should be noted that, of the sources encountered in the litera- 
ture survey, published data measurements found subsequent to those 
used herein agree to within experimental error. 

(2) Since some of the *'best” data available are generated by 
polynomial correlations of numerous sets of data (see belox^7 for 
Cp and k) , concern was experienced regarding the capability of a lov; 
degree polynomial to be able to predict accurately the behavior of 
higher order derivatives. For example, v;hile a second degree poly- 
nomial equation might adequately fit a set of data, the constant 
second derivative analytically derived from it may not be realistic. 
Erratic behavior in the derivatives can also be expected if the degree 
chosen is too large. Evidence of these two cases has already been 
presented in Figure B.l. Using the smoothed cubic spline as an 
approximate standard v;ith the data set of Osborne, ct.ai. [31] (judged 

to be one of the most reliable [30]), and knowing that only the 

dc 

first derivative, required for tliis analysis, one can 

see that the tliird degree polynomial of Touloukian, et.al.[30] 
adequately approximates the required behavior over the temperature 
range of interest. V;hon second derivatives are required as v/ell 
(such as is the case for p, p, and k) , the above j>roccdurc remains 
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basically the same, although more involved. For example, again 
using as an indicator the smoothed cubic spline v;ith data recommend- 
ed by Powell [32] (also judged to be among the most reliable avail- 
able for k [33]), first and second derivative behavior is found to 
be best approximated when thermal conductivity is defined by a third 
degree polynomial . 

In short, when an analysis of existing data uses polynomial 
correlations which w^ould obviously not predict behavior of higher 
derivatives needed herein, information characteristic of that ob- 
tained from the polynomial fit is used instead. 

(3) Regarding the required property variation of water, the 
following information is relevant: 

Specific heat , c^. An analysis is made by Touloukian, et.al [30] 
of twenty sets of data for the isobaric specific heat of liquid 
water. For the saturated liquid, recommended values are computed 
from 

Cp(cal g“br^) = 2.13974 - 9. 68137x1 0~\ + 2.68536xlO“^T^ 

- 2.42139x10“®!^ (T in "K) (B.l) 

for 273® < T < 410®K, This equation is found to fit tlic data 

analyzed with a mean deviation of .14% and a maximum deviation of 

1.83% (for one data set - all others differ by no more than .5%). 

Of particular interest for this study, reference [30] indicates that 

c as calculated from the above equation ’’should be substant.ially 
P 

correct to v;ithin one percent” bclov; 400®K. The ability of this 
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cubic to predict first derivative behavior has already been dis- 
cussed. 

D ensity, p . Four previous works are evaluated by Gildseth, et.al. 
[34] and compared with their experimental results. Recommended 
values for the density of air free water are generated for 
0® < T < 80®C by a modified form of the Tilton and Taylor equation 
(i.e., the inclusion of the exponential term which is negligible 
for T <40°C) 



p(g/ral) == 1 - 



(T - 3.9863) ( T + 288.9414) 
508929,2(1 + 68.12963) 



+ 0.011445 exp(- 
(T in ^C) 



374.3 

T 

(B.2) 



) 



to fit experimental data with a mean absolute deviation of 
.7x10 ^ g/ml. \duMe the authors do specify these values to be un- 
certain in the seventh decimal place for the lower temperatures, 
they do indicate that at least above 40^C uncertainty is in the 
sixth place and increases with temperature. Since this equation 
does reproduce the densities tabulated in the International Critical 
Tables [35] for the r^mge 80® < T < 100®C (to v;itliin the indicated 

uncertainty), it is extrapolated in this study to cover the full 
temperature range, 0® - 100®C. 

Dynam ic visc osiLy, p. The experimental measurements of ICorosi, 
et. al. [36] are compared with ten previous, similar sets of data. 
The authors propose a correlation 
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log y = ~ 1,64779 + 



262,37 
T - 133.98 



( y in cp and T in 



to fit their data over the temperature range 20®< T <150®C v;ith 
an average deviation of 0.17% and a inaximum deviation of 0.49% at 
40°C. However, since this data seems to be characteristic of 
the mean of the other compared data for a specified temperature, an 
alternate formulation is used. Interpolated data is generaged 
instead by using (for the same teiTtperature range) an equation 
recommended by R. E. Manning in reviewing the paper: 



This equation is found to represent Korosi’s e>rperimental data to 
v/ithin 0.05% average deviation. 

Thermal Conduct i vity , k . Of the more than sixty sets of experi- 
mental data available on the thermal conductivity of liquid v:ater, 
Touloukian, et.a], . [33] evaluate seven as being the most reliable. 
They then correlate these and generate recommended values for 
273.16^ < T < 4l3.16°Iv with the second degree polynomial 



This equation is found to fit the data considered with a mean 
deviation of 0.24% and a maximum deviation of 0.82%. As indicated 
in previous discussion, however, a second degree polynomial appears 




1. 37023 (T - 20) + 8.36xl0~^(T - 20) ^ 
109 + T 



(B.3) 



(T in °C and y (20°C) = 1.002 cp) 



10*^k (cesu) = - 1390.53 + 15.1937T - 0. 01903901’^ (T in '’K) 



/ 
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to be inadequate to describe the behavior of the second derivative 
in the temperature range considered. Thus, for application in this 
stability study, a third degree, least-square polynomial is fit to 
an earlier correlation of therraal conductivity data by Pov;ell [32] 
with the following resTilt 



k (inwat.tr> cirr^K - 9.901090 + .1001902T - 1.373892x10 



+ 1.039570xl0~\^ 



(T in '"K) 



(C.4) 



Values generated by this equation deviate from the preceding one by 
no more than .25%. 

Due to the nature of the formulation of the mean and dis- 
turbance equations, it was found convenient to calculate the pro- 
duct of the density with the viscosity and thermal conductivity, 
p |i and p k, rather than these properties individually (see equation 
2.20 for example). Thus temperature variations of p, p y, and 

p k as well as their teraperature derivatives required for this 
analysis are sho^^Ti in figures B.l through B.5. The ratio of the 
second to the first viscosity coefficients is computed using 
equation 3.3 of reference [37, p.47] and the data of Pinkerton [14] 
and is included for reference in figure B.6 (see section 2,1), 

Since the stability characteristics are found to be sensitive 
to assumed property variation v/ith temperature (see Chapter IV) , a 
comparison is made betv7een fluid properties calculated by equations 
(B.l - B.4) (those used for this j nvostigatl on) and by the 3 cast- 
square polynomial fits of Kaups and Smith [7] (used by V/azman et.al.) 
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Deviation of the latter from the former (of the form 

/OO j for example, v/here po is the reference quantity 
evaluated at ~T - 32^F) is Siho^vm in Figure B.7, 



kAUpS etol 
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Figure B.3 First and second tcnperaturc derivatives of p as a 
function ol tcnpcraturc . 



i 



(rld)p 



126 




Figure F.4 Fir^l. and second l:c>inpcral:urc derivatives of o\x as a 
function of Lenperature . 
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Figure 13.5 First nncl second temperature derivative of pk as a 
function of temperature. 



( pk) 
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Fl{>urc B.6 Second- to-f irfit viscom'ty ratio for water «as a function 
of teinperature . 
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Figure B.7 Deviation between the property-toiTipera turn variation of 
Kaups , et.ai [7] anci that specified in Appendix B.2. 



APPENDIX C 



MEAN FLOW SOLUTION VARIABLES 



For reference purposes, the Mean velocity, tenperature and 
pi'opcrty variation through tVie boundary layer, as well as their 
n derivatives coMputed b}’' the program are plotted in this appendix 
for the flows eonsidered in this investigation. 

To transform this information from the n coordinate system 
to a more physical one, Figure C,1 has been included. If, for 
example, the reference normal scale is denoted as y^ (i.e, y - , 



In this formulation the advantages and disadvantages of using n 
(rather than y) as the independent variable may clearly be seen. 

Advantages 

(1) Using the mean density variation, the physical step size 
is ’’automatically" adjusted \/ithin regions where the temperature is 
significantly different from its free stream value. This menus that 
although the n step size remains fixed throughout the entire in- 
tegration range, the actual distance traversed v;ill be varied 



o 



then 




or 
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constantly within the Ihermal boundary layer. The larger the wali- 
to-free stream teinpe,rature difference (and so, the more severe the 
therroal gradients near the V7all) , then the smaller will be the 
.initial steps from the wall. Admittedly, liquid v/ater*s density 
variation with temperature is not as significant as might be desired 
for such a formuJ.ation, but for large temperature differences, it 
is still sufficient to permit selection of a constant step size out- 
side the thermal boundary layer v;hile adequately describing varia- 
tion within it, 

(2) If information obtained from solution of tlie mean flow 
equations (2e20) is to be used directly in defini ng the variable 
' coefficients for the disturbance equations, then both must be 

expressed in terms of the same independent variable. If the latter 
is to be solved instead with y as the Independent variable (to 
compensate for the first disadvantage discussed below), then a 
constantly varying step 




must be used within the thermal boundary layer. It v/ould seem to be 
inconsistent (particularly for large temperature differences) to 
solve the mean flov; equations in terms of tlie Hov;arth-Doro tnitsyn 
variables (eqn. (2.1‘'0) ^nd then assume the mean density to be 
constant (p = 1) to solve the disturbance equations (as apparently 



mmm 
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Wazzan et.alo [ 3-5 ] have). Such an assuraption should raise 
questions concerning the accuracy of the variable coefficients for 
the latter, especially in the wall region where density changes are 
■greatest, 

Disadvantag es 

(1) Outside the thermal boundary layer, the n step size is 
tiraes that of the equivalent one. For the temperature 

range considered r, . , is alv/ays greater than unity (T < 213^‘F with 
“ 60°F, B = 0 for ) so even when the displacement thickness 

is chosen as the reference length (y^ == a larger n step would 

actually be required than for a comparable y integration, Ihis 
means to obtain the same accuracy, more steps would be required with 
X] as the independent variable, 

(2) An extra calculation is required to transfonn back 
from the q to y coordinate system. 

It was decided that the advantages outweighed the disadvan- 
tages and so the problem V7as formulated and solved using q as the 



independent variable. 







PlRurc C.l Rolatlonfihip bc-.tv/ccn nand for variou;; wall temperatures. 
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Figure C.2 Variation of witli changes in the wall tomiicratnrc. 
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Figure C.3 Second velocity derivative evaluated at the v*all as 
a function of the wall tcir.pcraturc . 
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Fip.ure C.4 First dprivativos of velocity and temperature 
evaluated at the uall as a function of wall tenperature, 
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Fip.ure C.5 Specific heat boundary layer profiles for various 
wall tcr.pcratures . 
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Fip,ure C,6 Density boundary layer profiles for various wall 
temperatures . 
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Figure C.7 BounHnry l^vcr profiles of ‘cho first derivative of 
density for various wall Lernperatures . 
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Figure C.8 F’onnclnip^ layer profiUs of the Jiacond derivative of 
density for various v/ali temperatures . 



141 




Fij'urc C,9 Boundary Wtyor profiles of tlie dcnsity-vl5c:o:^uty profluct 
for various v;aii tenpe raturos . 
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Flj^ure C.IO Boundary layer profile:? of the first derivative of the 
density-viscosity product for various v/all torperatures * 
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FijMire C.ll Bouriflary iayor profiles of the second derivative of the 
density-viscosity product for various wall tenpera ture.s , 
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Figure C.12 Dounrlary layer profiles of the dons Iry-thernial conductivity 
product for various v;all temperatures. 
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Figure C.13 nou.ndnry layer profiler; ol* the first <i(>rivativc of the 

density-therTPal conductivity product for various v;ali rer.porri Cures , 
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FifM.ire C.14 BoundiTv lay^^r prolljcr; of Lhe second derivative of 
the deiisi t y-tberrtal conductivity product for various waJi 
temperatures , 
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Figure C.15 Boundary layer velocity protilei; for various wall 
tenperaturcs . 
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Figure C.16 Boundary layer profiles of rhe first derivative of 
velocity for various wall tempevatures , 
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A 

Fipuro. C.L7 Boundary layor profiles ^or the second clGrivativo of 
vcioclty for v<irlous walJ Lenperatures • 
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Figure C,18 Boundary layer temperature profiles for various 
wall temperatures . 
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Fi^.urc C.19 Boundary layer prcCilc;^ of the first derivative of 
tCTTipcr/itiurG foT Vcirioijs wqII tcnncr^turcs . 
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Figure C.20 Boundary layer ]>rorileis of; thn, second derivative of 
f.cnperature for various v/n 1 1 tionpeiMtures . 



Appendix D 



The Computer PrcHxrain 
^ ^ Des cri ption 

The progxaia \nritten to solve the stability problem forxrulated 
in Chapter II consists of a mapping processor (BLSJIAP) , the main 
program (BLSTAE), and eleven subroutines, all of which are coded 
in FORTRAN V. A brief description of each, along with its role 
in the overall computational scheme and references to numerical 
procedures (where applicable)^ is given in Table D.l. 
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TABLE D.l Continued 
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D.2 In put 

A description of data required to successfully execute the 
program described previously is now given ii\ Table D.2. Further, 
the required foriviat for this data is illustrated by a sample data 
deck in Figure D.l. Note that data input is arranged so that the 
first card lists only the logical control variables for the 
program; the second, all mean-flov; related information; and, the 
remaining, disturbance flow and stability related information. 

If only the mean flow solution is sought, then only the first tv/o 
cards are required: all others will be ignored by the computer. 
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Program 




Card Field 


Variable 


Description 




Prop;raTn Control Pararaeterv^ 


1 


1-30 


NOTE (I) 


A series of logical variables identi- 
fied in BLSTAB which controls the 
operation and output of the program* 




Mean Flow 


Parameters 




2 


1 


TWALL 


V/all temperature, (specified in 

degrees Fahrenheit) . 


2 


2 


Tir^F 


Free stream temperature, (specified 

In degrees Fahrenheit) . 


2 


3 


BETA 


du 

X e 

Falkner-- Skan parameter, M = ^ 

e 


2 


A 


EPSILN 


Specified boundary layer edge para- 
meter, vrhere uL = 1 - EPSILN. 

0 


2 


5 


H 


Step size for mean flow equation 

calculations, h . It is half that used 
m 








for the disturbance equation integra- 
tion, h^, and must be greater than or 

equal to /700 where is the 

position defined as the boundary 
layer edge. 




Stability 


Parameters 




3 


1 


IlTJfAX 


Maximum number of iterations per- 








mitted to compute an eigenvalue. 


3 


2 


ANGLE 


/\n orthonormal basis is reformed if 
the angle between any two solution 
is degraded to less than this 
quantity . 



TABLE D.2 Description of Input Data for I'ropran BLSTAB. 
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Card 


Field 


Pirogram 

Vai'iable 


Description 




Stability 


Parameters 


(cont « ) 




3 


3 


EPSLK 


Indicates the dcp;ree of refinement 
of eigenvalue estimate* Iteration 
is terminated when |f(ITER )|5 








, RE (EV (ITER)) I 

' ~ RE(EV(ITi:R+i)) ' ’ 










1 Am\G(EV(ITER)) 










AIMAG(EV(ITER+i)) 










are all less than EPSLN. 




3 


4 


NEV 


Specifies which stability parameters 
are to be fixed and which are to be 
varied as the complex eigenvalue, 

EV* 








EV(ALFAR, OMEGA) 

2 EV(ALFAR, ALFAI) 

3 EV(ALFAR,RE) 

4 EV (OMEGA, RE) 

5 EV (OMEGA, ALFAI) 

6 EV(RE, ALFAI) 


A 


1 


IM 


Number of eigenvalues to be 
(maximum of 10) 


computed 


5-14 


1 


REDLS(I) 


Reynolds number, • 




5-14 


2 


ARDLS(I) 


Real part of the wave number, a* 6* 

K 


5-14 


3 


AIDLS(I) 


Imaginary part of the v;ave 


number , 








a '* ^6^' • 




5-14 


4 


FREO(I) 


Frequency , 










All of the four proceeding quantities 
are the best rstiiiiate of the Itli 
eigenvalue to be commuted. 



TABLE D.2 Continued 



164 



^ ^ Prog ra m Listing 

A complete listing of the FORTRAN V source deck used to 
generate the results in this investigation is given in this section. 
Note that v.^here possible 5 the routines have been written in general 
so that application to systems of order greater than six is 
possible imriiediately v?ith perhaps merely a change in the DU'IENSIOII 
stertements . This is particularly obvious in the routines 
associated with boundary layer integration > ADAMS and DIFFEQ, and 
V7ith the Gram- Schmidt Orthonoriaalization algorithm, GSORTH, 

EGNFCN, and ORTHCR. 
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29 
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33 
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c 

c 

c 



uOT£(29) 

..O*iX(30) 



M “ b 



PKliU AKbL£ fi^TVLEi: SOLfJ VlCTOPS ON FIRST iTERAflCN 

LObELL^b L^NAflOiJb. IR FALSE.. USE LOUATIGNS OF bAZZANf 
0<M>i UivA f A;.D SMIfH. 



N = 3 



K..;f,X = Vul 
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Alpha = (beta + i.o)/ 2 .o 

R = Tli- l /'lOl .69 
S = (T..Al,L-n;jF)/49.l.69 



U 



C2 = (T...M.L-TIt.K).f5./9. 

C.ALL TLrvAK(O.OfO) 

IF (^.0Tf:(9)) .VRITF (6»1(:9) 

11. -EPS U-rl 

} FOKSAT (' STEP SIZE ='riU6.a/' falkher-skah parameter 
1' V.ALL TEimPCR/MURB -'rFj.0.4/< FREE STREAM TEMFERATUR.E 



! U BET A - rv'ALL-459 . 69 f T 1 tlF~(|59 . 69 » PR 1 NF f 



»F10.4/ 
:» rF10.'4/ 



couriuARY l-ayer edge UF.TERMIHE 

1 



Id 



2» FREE S1REAM PKAMOTL MO. ='fF10.b/‘ 

30 by U =' ?R (1.7/) 

IF (KOfEdi) .AMO, .MOT. MOTE(D) (iO TO 

Call i so. i ( n » eps i lm » fv all . 1 1 mf r etadls ) 

ETADEL = bLPAR( I rKTOT) 

IF (MOTEC.;)) VKITE (f.jlOl) ( K r ( BLP ',R ( I » K ) » 1 = 1 f 7 ) r K = 1 r KTOT ) 

FoliMAT (/:. i;<r ' * • ,15Xd , , ' .37X* * < ' d6X» ' , » »/5X; «K‘ f 1?X» 'ETA* r l6Xr 
1 f 17 a» , ISXf 'v.' fl7X» ' TriETA* » 1 2X d THETA '» i2X d THETA V 

? (.(t)»7C.10.B) ) 



IF (;nTL('.)) -.rite (6,102) KTOT.ETADEL 
io2 Format (/' i-tot =-.i 4/' ETAfoa. =',ei3.o/) 

If (i.OfEd)) V.RiTr (6,lo6) ETADLS 
Ijo (/! cTA SOcSCRli'T OLLTA STAR =',El6.S/) 

v.riTi: !b) r'M.Tiir ,ET>-r LS,LT aMEL-KTOT ,M0TE (24) ,n,T,-.'ALL»TliJF,BrTAf 
IE, SILM' ( (.'H AA..(J ,K) rX = l,7) .X-ltKl'Ml) , ( (CU'K) , I = ir 14) ,K = 1,KT0T) 
lA (..'Gi';n)5 SToi' 

00 r 0 2 

1 RE«U (u,FMO=cT) PRIMF, etadls, ETA 0a.,KT01 , H20VAR , HA , TVMLLA , T1 1 IFA > BET 
1A,.,EFSA, ( ( ;LPAR( C-r.) ,I = j ,7) ,K = 1,,. lOr) , ( (C( I,K) ,1 = 1,14) ,K=.l,K ror) 

IF (iiOi .</') .Ai/J. H23VAR) GO TO i'J 
Ir (••:: !K (2 : ) .Of . H207,.R) go to 11 

10 IF ( iBlh! ( (^-'.-n) ' 1 .CL6) .£0. 0 .•■..•-•j. TMT(((T VAl LA-TWALL)-!()'Ii,;-A" 

iTir.F) 1 •.••iiLrA) + (i-:"SM“EpsiLi;) ).) .EG) .eo. u) go to 9 

11 ,i ;ITE (l, ;.1,-) 

il7 r\,'i..-.',T ' 4- < ) /O/ . • fHE RCCUF,;. rro riFORMATlOM IS mot the SAl'E 

IAl Tm.-.I Sl.il ED IM Uf.n t . '/ix,0b( > , • )/) 

G T Op 

9 1- (:.-.j:E(.-,, ) (RJT,-. (6,lJi) (K, (I'Ll ;.i(I,r) ,1 = 1,7) ,K = l,KTOr) 

ip (.j'/rp C. ) ) ( 4 - ’.,■ 2 ) KIOT,'. l/l'•i•:l. 

If ( ■/., I .. ( / ) ) ..'<11t-, (■,!,<-') liT/ BLS 

2 ( ., • f. A 1 ) I I A; .'Aa , 4. , i.f^SLM , f-.t-V 



••Of P \1 
If (.J.ji 



' t , i* I (1 . I , t„ ) , I 1 1, i 

_')) v.Rdc. (K, (C( 1 ,!• ? , 1-1, 14) ,K.-;l ,KT0T,5) 
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loa 

ruy 

i.VQ 
111 
112 
11^ 
li^^ 
Ub 
11 b 

IIV 

lib 

119 

lao 

lai 

laa 

lab 

la^ 

12b 

126 

i27 

126 

129 

130 

131 

132 

133 

134 

1 35 

136 

137 
130 

139 

140 

141 

142 

143 

144 
14b 

146 

147 
140 

149 

150 
Ibl 

152 

153 

154 

155 

156 

157 
15d 
159 
lUO 
161 



lu3 FOhP.AT ( ( 14/ 14ti0*3) ) 

IF (i;urE(3U)) WKIJE (6/107) 

IF (.;;0T. .^GTL(oO)) \-.R ITc ( 6 / 100 ) 

IF (*A'OT. n01L(^3)) .ViUTC(6/U0) 
lj7 F01-^5^T (/' KE5ULTS FO.< FFUATIOUS OF LO.vFLL 

lOO FoKMAT (/f F;£5ULTS FOf^ FQUATIO:jS of WAZZANr OKAN'.URA/ AMO SMITH 

1 

110 F..AMAT (» V.ITHOUT IFMFORATUFE FLUCTUATIONS 

IF (N0T£(6)) WKITF (6/104) ATiCLF ? I TRMAX / FPSLM 
iu4 FoKMAT (/» MIN At, OLE ALLO..LO «'>FTV;l'N SOLM VECTORS =’fFlO,4/ 

I* i'^AX i^C* iTEKATlCi^S PERMITTFO =»rl3/ 

2* CONvLPGlNCE CIUTLRION For SIGEfJVALUES ::'/El0.4) 

H.:,7EP = -2.0=<4I 
K:.:AX2 - (Kr:AX/2)+l 

c EbriM.-FfL: i?E/0’;^0/w AMU ALFA TO USE IN E«V, CALCULATION 

READ (5,203) FiF/ (HtCLS(I) / ARDLS ( D /AIDLSd) ?FRE0(1) d^l/MM) 

2j3 FoF<MAT ( 14/ (4L16,0) ) 

Do 3 I“l/KiiM 
NjTE(2l) = •TRUE, 

KL = RLDLS(I) /ETADLS 
Auf-Al = AI::LS( I )/ETADLS 
0..-.LGA = FREO(I)^RE 
ALKAR “ AKoLS( I) /ETADLS 



IF 


CCV 


• EG. 


1) 


£V(1) 


i: 


C; iPLX( ALFARf0f.;e<3A> 


IF 


(iJc-V 


• uO • 


2) 


EV( 1) 


:: 


CMPLX<mLFAR»ALFAI) 


IF 


<:;f.V 


• EO, 


3) 


EV 1 1 ) 


= 


Ci iPLX( ALFAKfUE) 


IF 


C;lV 


• EG. 


4) 


£V ( 1 ) 


:r 


C!-,Pl.X(Or-.tlGAfKF) 


IF 


(,\'CV 


• EQ • 


5) 


EV(1) 


n 


CP.PLX ( OMEGA f ALFAl ) 


IF 


(NCV 


♦ EO • 


6) 


EV ( 1 ) 


- 


Ci-',PLX(P.E»ALFAI) 



Ev(2) ~ (1.0i/0,0)-i^r.V(l) 

Ev/(3) ~ CMPLXn .0f-5^REAL(CVU ) ) p1,02^AI !AG(EV(1) ) ) 

Call ulSTEGULV/ iTfAMAX/hSTEP/ A ij 6LF/Y/EPSLN/ PR I MFNMEV^$3/Pr 2D I ST ) 

FREO(i) :: OiMEOA/RF 

ARuLS(I) :: ALF AR^i. T ADLS 

Ri:OLS(l) z i^EF^bTAr.LS 

AIDLSd) = ALTAI I'M alls 

C,< = 0!-I|:'}AvALF/‘R/(A!.F.U. ^=?^2-^ALFA[^’^2) 

C / — *~0; lGA'^ A f / ( AT Far ^^2 ' ALF i \ I ^'2 ) 

IF (uOTl(o)) (ofJoG) AR0L5CI) iAIDLS(I) rkEOLS( 1) /F};E0(I) 

l^CR/Cl 

lu5 F^RfM/r (/» ALFAR =:»/FlucS/» ALTAI ='/Fl6*3/’ RE =»/E19,A/» FRLO n: 
Idclo.o/* CR =^L.l9.F/* Ci -»*E19,0) 

IF (••.v)(E(2) .OR» .:U0T. ‘iOTr(2l)) 00 TO 3 
K.'TOT XT01/2M. 

IF C.:.UT. r.0TF(l5) •CR, .NOT. NOTl(30).) GO TO 5 
.v.% ; T r: ( o / 1 j, 2 ) 

lj.2 {//» KP» r'-X/ N'luNSn Y hLNCTUATiO:/ //:</ ♦VISCOSITY FLUCTUATio 

IN TPLUPAL COi.UUCTlVITY FlUCTUATIG:^ ASSdLNFLU) AI>b(VISFLU) A 
EdSdCrLU) ' ) 

0;> ^-'i/KPToV 

TiF:P “ FL ^.u;(5/R^:.'‘l) 

ChlL rE‘:JA;>(TL? P» 1 ) 

Ii- (.,oN. (,'.,>) ) KLUcro(R/.d r. 0 ( D ♦'GP ( 2 ) ^ Y ( 3 / K ) 

IF (.UNT, NOTLCFo)) rLUClG(2/K) - (O.O/O.C) 
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iba 

16Ji 

Xo‘f 

165 

166 
167 
lob 

169 

170 

171 
1 U 
176 
17n 

175 

176 

177 
17b 
179 
liiO 
161 

16 5 
164 
13b 
186 
187 

138 

139 
190 

un 

192 
196 
I9't 
195 
] 9o 
19/ 
1 98 
j 9 ‘ 
auo 
aoi 
aoa 
2 1? 6 
2;}4 
29 5 
auo 
;:a7 
a 08 



FbUCTOti^K) G( 3)’/GP(3) ^Y(3 m\)/G(2) 

FuG0TG(4i.s} •* G(4)*bP(4)^Y( 5^K)/G(a) 

Ir (lul (T..;\LL--TIr3’) 0) V.'iUTE (6^11) K ^ ( FLUCTG ( J r K ) ^ ; 4 ) > 

i(LA;iS(FLUCi'3(J/A) ) rjr.a,4) 
ill fOh’v.M ( l4ro{fl4,d/2X^3E.14,8) 

4 C(j,nTI3LJF 

5 IP (pi JUT. .vOTF. (16)) GO TO 3 

C TFk.jS of LiNF:-,0Y B;\LA.iCF EQUATION FOR NEUIKaL D 1 5TURG ANCES » 

C FUN0(i,K) aEY.JOLOS STi^cSS PRODUCI lOU 
C FUNv,(3rK) ~ r:..3LUCTl0r-. OIT: fO V15C0SITY FUUC TU AT lOf-iE 
c fl;)o(o,k) = (\-;>::b5URc ee'Ergy pROoucTiori 
C. Fu.'-;C(4,k) r: j 1 35 1 PAT ION 

C Fu.'lc: (5;K) :: jlbSlPATlOii Dl^E TO r.OrJCODS TANT MEAN Af^'D FLUCTUATION DEN, 
C FUNc(GrK) rlvOUUCTlON OLL TO FEAN VISCOSITY GRAOIEl.T 
«v.uTE icfllD 

ii4 FOH?wU (/» ENERGY PA'OOUCT lorj AND DISSIPATION TERMS FOE NEUTRAL STa 
lUiLXTY (ALFAl::inO) ♦/» KP RE STRESS PROD VISC FLUCT PROD PRES 
as ppoD’ roAr ’DissipATiorN »o » »dem dissip VISC gkao prod*/) 

Do 6 K^l^KPTOT 
L - 2tK-l 

FoKC ( 1 mo :: -ULPAR ( 3 . L ) Al.F AR + RE-^^RvEAL ( Y ( 2 > K ) ^^COf ! JG ( Y ( 1 FK ) ) ) /C ( 1 > L ) 

Xp (pijOTe :.oTE( 50)) FlJIX(l.'K) = FO. IC ( 1 F<) ^ C ( 1 f L ) 

FO:^:C i^fK) - -DLPAR ( 3 ? L ) 4 ( i^E AL ( F LUC TG ( 3 > K ) ^CoN JG ( Y ( 5 f K ) ) ) /C ( 1 a. ) 

1 2.4- A }. I 'AG ( Col>i'UG ( Y ( 1 M\ ) ) ^FLUCl G ( 3 f K ) ) ) 

IF (,i\‘OTp ,.i01c(30)) FUi;C(aFr<) - 0«0 

FunC (3»K) - -ALFAF4^r<E^FLAL(COi‘!JG( Y(4f K) )=«■'(( 0 1 . ) ^^FLUC TG ( 2 ;K) ^ 

i (f.LPAR'2^u)"-a:)4C ( I rL) 4C(14f L)/2, t.y ( 1-K) ) ) 

if (,f;Or* i.’orECSO)) r-uE.rfSM;) = o*o 

FuNC ( 4 » K ) ^ -c ( 2 f 1. ) :^ ( AI.FAR^^^ MT + Y ( i MO -J CO’ 1 JG ( Y ( O K ) ) ^: C ( I f L ) 

1C. JLJO ( T ( 5 F iO ) ( 5 M<) /C (1 ^ L ) ^^^2 + 2 * -MM-FAR^ + e^-A IMAG ( CON JG ( Y ( S MO ) * 

2Y(i fK) ) ) 

Fu;:c (Sf K) - a . (aa.) c (7 m.) ^'FLi^aR (g.»l) / c ( 14 pl) ■^=ai.far + 4 ^ 4 . ( ( ( 

1G..PA’; (F »L) -CR)^C ( 1 pL) ) r-F LUCTG ( .< M ) ^- COl -'JG ( FLUCTG ( 2 r ^ ) ) V (C (14 Mj / 
22. )^ (24r ( I f lO 4-C0A'jG( Y ( 1 / K) ) + ( DLPArO 2 r L ) "CR ) ( 1 4 M. ) 4 C ( 1 f L) 

3 A I nAG ( Lonjo ( FLOG I •' ( 2 M< ) ) 4 r ( 1 r K ) ) ) 

U (m: 0 i\ .,oTl:(30)) fu:c(SwK) ~ o*o 

Foj-.C (oMO Al.FAl- ^=^2 NT ( 2 r L ) /C (1 f L) -t ( C ( 6 M^ ) VC (14 ? L ) /2 O 'VAli-lAGC 

Y ( 1 ME ) ) v Y ( 2 r K) ) 

IF (CR ,GT. LLPA7(2M0) GO TO 7 

IF (Ci\*GT o4 P.\r{'dtL-2) ) .UTE(6r lib) K-1 p K f BLfMAR ( 1 ) MTLPAfUl^L) 

i^;> Foi-r..\r (» vriE critical point is fere wmeRE' » ilf *<k<» f io» » or 

lfO*cf*CET.X’M-4.2) 

/ OIL A; (njNC(J/K) f J = 1 f6) 

1)3 F,,Li AT ( I4f tiEUj.8) 

6 Continue 
3 Co.iTINUE 

ISO 



1 

b 

b 

7 

a 

9 

10 

11 

la 

li 

14 

lb 

16 

17 

18 

19 

20 

21 

22 

2^ 

24 

2b 

26 

27 

28 

29 

bO 

31 

32 

33 

34 

3b 

36 

37 

36 

39 

40 

41 

42 

43 

44 

4b 

46 

47 

4b 

49 

bO 

51 

52 

53 
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SULROUl INE MFEQN ( H r EPS I LN r TWALL f T I NF r ET ADLS ) 

DlMENSlOfJ r (15) fDY(ib) fETAEND(4) r TEST (4) r 8 ( 7 r 70 1 ) r C ( 14 / 701) 

1 f GN ( 4 ) ^ Gr^N ( 4 ) 

Double precision YfOYfhfETAfXf ZrOEEXfOELZ^SUMODOr sumevn 

LOGICAL note 

COMr.ON /8LKl/ALPHAfBETArPRlNF/BLK3/G(4) f GP ( 4 ) r GPP ( 4 ) /BL.K4/B f C/ 
lDLKb/N»MfKTOTfK(MAXfKMAX2/DLK8/NOTE(30) 
external DIFF 

DATA ND/2/f (ETAEND(I) r 1=1 > 2 ) /2 . r lO • / f (TEST( I ) f I=lr2)/1. f l.E-8/ 
1 iETEST/1*0E-12/ 

X = 1*UD0 
2 = -1*0D0 
1 = 1 
Mark = i 
K = 1 

3 ETA = 0*0D0 
Y(l) = 0*0D0 
Y(2) = 0*0D0 
Y(3) = X 
Y(4) = l.ODO 
Y(b) = 2 
Y(6) = 0*0D0 
Y(7) = O.ODO 
Y(8) = l.ODO 
Y(9) = O.ODO 
Y(10) = O.ODO 
Y(ll) = O.ODO 
Y(12)= O.ODO 
Y(13) = O.ODO 
Y(14) = O.ODO 
Y(lb) = l.ODO 

call AUAMS ( 15 f H r ETA r 0 f Y f D Y r 1 r D IFF ) 

GO TO (4f30) f MARK 



30 B(1»K) 




ETA 


B(2»K) 


“ 


Y(2) 


B(3*K) 


z 


Y(3) 


B(<+*K) 


z 


UY(3) 


0(5»K) 


z 


Y(4) 


ri(6*K) 


z 


Y(5) 


B(7»K) 


z 


DY(5) 



C COMPUTE MEAN PROPERTY VARIATION THROUGH THE BOUNDARY LAYER 
Du 1 I=2f4 
Gn( 1) = GP(I)*B(5fK) 

1 GlJfid) = GPP(I)*(B(6fK)t*2) + GP(l)*B(7fK) 

IF CriOlE(13)) WRITE (6d01) K r G (1 ) f ( G ( I ) > GN ( I ) ( I ) ► GKN ( I ) *G ( 1 ) r 1 = 

12. q) 

iol Format ( 14 . 10 E 12 . 7 ) 

C COMPUTE MEAfJ FLOW COEFFICIENTS FOR DISTURBANCE EQUATIONS 
C FOR MECh PHESSURErSET VIS2=-2/3. FOR THERMO PRESSURE^ SET VlS2 = 2.1 
VIS2 = -2.Z3. 

DVIS2 = 0.0 
Ca.K) = 1.0/G(2) 

C(2»K) = G(3)/G(2) 



b4 

bb 

bb 

57 

bd 

by 

60 

61 

62 

6b 

64 

6b 

6b 

67 

6d 

69 

70 

71 

72 

7b 

74 

7b 

76 

77 

7d 

79 

60 

81 

82 

8b 

84 

8b 

86 

87 

88 

89 

90 

91 

92 

9b 

94 

9b 

96 

97 

98 

99 

iOO 

101 

102 

10b 

104 

10b 

10b 

107 
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C(bfKJ 

C(4»K) 

C (bf K) 

C (6fK) 
1DVIS2 
C(7»K) 
C(8fK) 
C(9»K) 
CClOfK J 
C(11»K) 
C(12»K) 
C(lbfK) 



= (2.+VIS2) ♦ (2.*GN(2)**2-GNN(2)-GN(2) ♦GN(b) ) -GN ( 2 ) *DVIS2 
= GN(2)*(2.+VIS2)-GN(3) 

= -GP(3) 

= (2.+VIS2)*(GP(2)*(3**GN(2)-GN(3) ) -GPP ( 2 ) ( 6r K ) )-GP(2)^ 

= -(2,+VIS2)*GP(2) 

= -GN(3) 

= -GPP (3)*B(6f K)*B(3f K)-GP(3) ♦BCUrK) 

- -GP(2)*(1.+VIS2) 

= PHINF*G(1)/G(4) 

- -GNN(4) 

- -2*0*GN(4) 



C(14»K) = 2**Gn(2) 

IF ( (l.-Y(2) ) .LE.EP5ILN •AND. Y (4 ) .LE.EPSILN) GO TO 13 

8 K = K+1 

4 call ADAMSdbrHrETAf IfYrDYr IrDIFF) 
b GO TO (6» bo) » MARK 

6 IF(ETA.LT,ETAEND(I) ) GO TO 4 

7 Bll = r(7)**2 + Y(9)**2 + + Y(10)**2 

Bl2 = Y(7)*Y(12) + Y(9)*Y(14) + Y(8)*Y(13) + Y(10)*Y(15) 

B21 - bl2 

B22 = r(12)*=^2 + Y(14)*«2 + Y(13)**2 + Y(15)**2 

Cl = -( (Y(2)-1.0D0)*Y(7) + Y(4)*Y(9) + Y(3)*Y(8) + Y(b)*Y(10)) 

C2 = -( (Y(2)-1.0D0)*Y(12) + Y(4)*Y(14) 4 Y(3)*Y(13) + Y(5)*Y(15)) 
DEM = B11+B22 - D214B12 
DELX = (C14B22 - C24B12)/DEM 
DELZ = (C2*B11 - C1 + B2D/DEM 

IF (NOTE(IO)) VyRITE(6rl04) ETAEND ( I ) r X ^ DELX f Zf DELZ 
104 FORMAT (FI 0 . 6 r 4D20 • 8 ) 

E = (Y(2)-1*0D0)*42 + Y(4)4*2 + Y(3)**2 + Y(5)*42 
X = X 4 DELX 
Z = Z 4 delz 

9 IF(ABS(DELX/X) .GT.TESTd) •0R,ABS(DELZ/Z) .GT.TEST(I) ) 60 TO 3 

10 lF(E,Lr.ETEST) GO TO 12 

11 IF (I .EQ. ND) STOP 

1 = 14 1 

Go To 3 



12 mark = MARK 4 1 
YEND = Yd) 

ETASTP = ETA 

IF (NOlEvlO)) V.RITE (6d06) ETESTrEfXfZ 
106 FORMAT (/» MAXIMUM ALLOWABLE ERROR IN ASYMPTOTIC SOLUTION =dEl3*6 
l/> ERROR IN FINAL ASYMPTOTIC SOLUT ION. = d E 13 . 8/ » FINAL ESTIMATE OF 
2: wall VELOCITY GRADIENT = d D16 • 8/22X p » WALL TEMPERATURE GRADIENT 

3 -» pDlb,8/) 



IF (N01E(13)) WRITE (6rl07) 

107 F0RfMAT(/37X» d d IIXp d » ' r23X> ddllXpdp’ p22Xp ' r ' dOXp * p p V2X> 'K» p 
13X» 'SPbC. HT. dbXP »UE(JSITY» p4Xp ’DENSITY’ pSXp ’DENSITY’ pSXp ' (DEN^VIS) 
2’ p2Xp ’ (DEN^VIS) ’ p3Xp ’ (LEN^VIS) dbXp ’ (DEN* TO d3Xf ’ (DEN*TC) dbX> 

3’ (DEN»rC) ’ ) 

GO TO b 

C ENSURE KTOl IS ODD FOR STABILITY EON K^iTEGRATIOfsl A(40 ETADLS SOLUTION 
13 IF •EG. K) GO TO 6 

KTOT = K 
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lua 

109 

110 

111 

112 
11^^ 

114 

115 

116 
117 
ilB 

119 

120 
121 
122 
126 

124 

125 
125 
127 
125 
129 
150 
161 
162 
163 

134 

135 

136 

137 
165 



C C;..LCU^ATL’ SUOSCJ^II’T OCLTA STAR = [:iADLS USiriO COiVPOSlIt: SIMPSONS 

If (i;;T (T.wM,L-TiNF) .MF.o 0) CO TO 14 

Dc.I«tSi = r.lA 

ETA'JLS - CVASrp - YENO 

GO ro 17 

14 5'JFOOD - 0,0 
sjf.Evn = 0.0 
X.^2 r KTOT-2 
I.ii = KTor-i 

IS K::3,IM2p2 

15 Sj.mOI'L ^ SUP.OIjU 4 cdno 

D 1 6 — 2 f I rt 1 p 2 

16 SiJiTlvM - SU’*-EVM 4 C(1#K) 

Ut.I-lNT = (0 ( i ^ i ) ^C Cl t KTCT) j-a, ■<^SUVOaD44.4SUMEVii) 4ETA/3./IMI 
eVAOLS - •jEf-lMT-YKhO^ETASTP-ETA 

17 Ir’ (MGlE(7)) SfOlTL (6rl02) 1 . -EPS 1 LN ? DENli 11 /ETADLS 

lu2 FoKf'iAT (/» (LOUNDAKY LAYEK TMiCKiiESS ( d FQ r 7 r » ) /D I SPL ACEM£| jT T}i 
llCKfT;:S5) =<?FJ2,a/) 

IF (.001, NOTE (12)) RETURN 
K|0Tn2 :: KIOT - 2 
Yi,o;.o:i ^- 0.0 

WfUTS (6d.06) YhOflOi-WBC 1 » 1 ) 

lu6 Fi/Rr-Ar (• Y/(01SPLACEMENT TH I CKNESS ) » r 6X d ETA » /2E20 . 8 ) 

DO 18 K = 1 ri\lgi::2»2 

YjmONOM = Yf •ONrjf-.+ ( 3 ( 1 r K-l 2 ) ( 1 r K ) ) ( C ( 1 M<) 4 4 , *4C ( 1 » K4 1 ) 4C ( 1 f K42 ) ) /6 « 

1/lTAOLS 

16 V/KITE (6r.l05) YTiO; iDM » B ( 1 # K+2 ) 

1C5 FuEMaT (2E20.8) 

REiORN 

EfcD 
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Uli'F 



* 1 


C 


S j jf\OUT 


k 




SUbKOJT 


5 




OvUtiLf: 1 






0i2.L.\bi' 


b 






6 




’ 


7 




CmLL tl: 


b 




00 1 X"- 


9 




Ga(I) - 


10 




1 G/.(I) 


11 




Do k 


1^ 




G J ( I ) = 


13 




G, *X Cl) : 


14 




2 GjZ( I ) : 


lb 




Orel) = 


16 




Dr(2) ~ 


17 




DY(3) ^ 


10 




Dr (4) = 


19 




Df(b) = 


kO 




Dr ( 6 ) = 


21 




DY(7) = 


22 




Di'(a) = 


23 






24 




DY(9) - 


2b 




OV(lO) : 


26 




KGX(I)'^ 


27 




Di (11) : 


26 




DY(IO) ; 


29 




on 13) : 


30 




l+2.0tY(: 


31 




DY(i4) : 


32 




OY(Jb) : 


33 




1 (G2(l ) 


34 




KrJUKN 


3b 




End 



i;(2) »Gi;X(2) ?gn?.(2) 
) »GP(4) f GPP(4) 



10 CaLCUI.ATc The. UlFFtnl-r.'iTIAU E0UMI0N5 
ClFF(i;TA>Y >DY) 
t: 1 A » Y . OY 

0‘i Y t lb ) » OY ( 1 5 ) f GX ( 4 ) ! 02 ( '; ) » 

/.4.M / ALPMA / pcr/v » p;< I I'iF/iU 
Y ( , ) 

••■•Vm.M rEFp» 1) 

1 f 4 

OPlDt^YCy) 

0PU)«Y(1<J) 

b I <1 

OP ; n * r ( 5 ) 

" CPP ( I ) % Y ( b ) 'Y ( y ) +6P ( I ) ty ( 1 0 ) 

= oPP ( I ) V Y ( b ) 4 Y ( 14 ) +GP ( I ) : Y (15) 

Y(2) 

Y (5) 

-oi i ( 3 ) +Y ( 3 ) - ! PFTA* ( 1 . 0/G ( 2 ) -Y ( 2 ) * 42 ) Y aLPHA+Y ( 1 ) 4Y ( 3) ) /G ( 3) 
Y(d) 

-O;. ( 4 ) *Y ( 5 ) "PHIP'F+G ( 1 ) *ALPHA*Y ( 1 ) yY ( 5 ) /G ( 4 ) 

Y(7) 

Y (»)) 

-^X ( 3 ) +b Y ( 3 ) ~GI I/. ( 3 ) 4 Y ( 3 ) -a; ( 3 ) * Y ( 8 ) Y ( bFT A Y ( 6X ( 2 ) /6 ( 2 ) 

2) 4Y (7) ) ~ALPHA+( Y(6)YY(3) i Y U ) 4Y ( 8 ) ) ) /G ( 3 ) 

Y ( 1 0 ) 

n -OX ! 4 ) *0Y ( 5 ) -GMX ( 4 ) 4Y ( 5 ) “Gl I ( 4 ) *Y ( J 0 ) “PRlNF*ALPHA-4G ( 1 ) * 

7 ( 1 ) 4 Y ( 5 ) Y Y ( b ) *Y ( 5 ) Y Y ( 1 ) 4 Y ( 1 0 ) ) /G 1 4 ) 

= MU) 

= Y(13) 

= ~C-2(3)yDi 
2) rY (12) )“.' 

1(1.5) 

= -o2 ( 4 ) +DY ( 5 ) -&:-;2 ( 4 ) + Y ( 5 ) -0! 1 ( 4 ) * Y ( 15 ) -PR I liFY'ALPI (A + G ( 1 ) * 



lI'P.A>MY(U) 



( 3 ) -GN ( 3 ) Y. Y ( 1 3 ) Y ( G FT A4 ( GZ ( 2 ) /G ( 2 ) 
Y(3) t Y(1 )yY(13) ) )/G(3) 



1 

• 2 

4 

b 

6 

7 

6 

9 

J.0 

11 

12 

lb 

14 

lb 

16 

17 

18 

19 

2U 

21 

22 

2b 

2H 

2b 

26 

27 

2d 

29 

3U 

31 

32 

33 

34 

3b 

36 

37 

36 

39 

40 

41 

42 

4 3 

44 

45 

4b 

47 

46 

49 

bU 

bl 

b2 

bo 
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s joivO^T It Jr. rE:;VAR ( i ^ v-ioc \ ) 

Lt^GlCAL ;.OTM 

C^k..NiO,: /•^L^2/•rOS'3x^F(4; ?C2/BUK3/G (4) ^ GP ( 4 ) f GPP ( 4 ) /BLK6/r JOT F ( 30 ) 

I /::c-K I /ALPi lA f bt.T A > P:n I PF 
DlKLPblO; m( 21) f.3(19) 

D„1A (M(t.) , L-i » cii )/l,'ir33oa9,0.f'07c:501»T..<?.a060^»0. 00592a >0.4C 159 j 
11,0. 23 o') 7 , 0 . 023'. Oo9 , '.5 = 259ovo0 r 137 . 1 9-530o3 >142.7970 /4 > b9 . 019762 

23,9,9a---U,--!,-2.5a201-.'6,6.7o062-n >3.734'1471>0.61A3'H7.73.3V690&.. 200. 
j7 ‘7455c, 197. < o04o ' 6b • c'6261<>t> , 7 . '1 7794 5ri/ 

D„ T 4 ( 6 I L ) , L- 1 , 1 : ) / . 1 3 r 7‘t > -9 , o 1 37E-3 > 2 . o6 J 3oE -5 , -2 . 4 2 1 392~8 , 

13 , 9663 * 2»u . 94 14 > 5.)(, 129 » 2 , o4 . 1 2963 , 0 . 0 1 1 445 , 3 i"l . 3 > 1 . 002 , 1 . 5 /023 > 
20,3bL-4 , 1 09. ' .43429 <43, -9.90109, . 10 019 32 >“1 . 073«92E"'( > 1 .n39570"7/ 

c ir- iN.,F.x = 0 00 I'uTh'S roc?: gtk£a:i vALues, Gir;i-' (kaupo s s;.;iTH form 

C COMPvTF-S OI.-.F/GFFF KOF uUAn’ilTTFs EVAOUA TLD AT T - 32 D 

C IF IiJ,,c-X > 0 CO'-FUTF f-!;('rF.KTY VARIATIOl;, G(I)> AT SPECIFIED TEMP 
C GU) = SPECIFIC MEAT aT COFSTAFT PRESSURE 
C G(2) u CEi;SlTY 
C 0(3) = GE..S1 V/tVISCOSITY 

c G(4) - r i-*( T hermal cofcuc (ivity) 

c ALL These cua,.tities are ucldim, l.r.t. reference values at riiiF 

T = S-!=Tl:-P4R 

TC = ( l-l)i-E73.16111 

TK = TC + 273, 1& 

If (TC .Lf. 1.0) TEXP - 0.0 

IF (TC .G£. 1,0) TfcX(’ = EXP(-B(10)/TC)/TC*-«'4 

IF (.liUT, mote (24)) GO TO 1 

0(1) = r ( 1 ) +8(2) <'TK+lt(3) ■7TA**2 t^B(4 )>->TK**3 

G(2) = l.“(TC“U(5) ) » + 2+ ( TC+9(t.) ) /b ( 7 ) / ( TC + B ( 0 ) ) +B ( 9 ) ‘HEXP >TC+ »4 
0(3) = 0(i) t3(ll)/10. + +( (3(12) +(TC-2,0. )+B(13):!‘(TC-20, )1‘-!<2)/(TC + 
lb (14 ) ) ) 

G(4) = G(2) + (B(16)+B(17) *TK+B(18)+TK‘»*2+B(19)+TK*i‘3) 

G\^ To 6 

1 G(l) :: A( i)-A(2)+T + / (3)+=T+*2 

G1.3) = A(4)+A(5)-»T-A(&)*T-#»2+A(7) >T + >i‘3 

0(3) = 1 . c/(.\(c )-A (9) ' r + A( 10)‘>‘T’+2-A( ll)4‘T + =i'3+A(12)*T*+4) 

0(4) = ",v(13)+A(14)-'1"A(15) ‘V + + 2I A(16)+T=+*3 

3 If (li.ufX .GT. 0) GO TO 4 
Do 2 1=1,4 

2 Glf.Fd) r- G(I) 

If (.•iClE(.j ! ) 1 P.UlF " 41 .04 + G(1)*G(3)/G(4) 

IF (,i:'-’f. ,.GTL(24)) pri, F = l3,66/(A(17)"A(lC) + r + A(19)*T**2~ 
ia(2'j ) .- 1 ‘fj 3 + ..(21 ) 4 ) 

RETU,-;. 

C C.LfU.E (7G( ( ) /ijH) /G( I ) AfiO ( U20 ( 1 ) /UM2 ) /G ( I ) '/.'MEfO; U = (T-TlfiF)/ 

C ( !..AL_-( ; ■-) 

4 if (.;.oi. :,.;irE(;'i) ) oo to o 

OF ( 1 ) = ( ,, ( 2 ) ■• 2 . u". ; P: '. i 3 . +B ( •■;)<■;». i +2 1 + C2/G ( 1 ) 

CF(2) = (-' TC-r ( L- ) ; : (3. TC-.2. F • ( 6 ! -’• ( 5 ) ) /B ( V ) ) 1 . -0 ( 2 ) +0 ( -^ ) +TFXP* 

1 1 c ( TC . 2> . 'I" ; ! 1 : • G2. • ( 2 ) / ( tc+h ( r- ) ) 

G••I•(F) :• ( ( (-U. “ ‘ I • . L-. ) ••2 . « • I'j) ) /'- (7 )-. •j(M) >!- ( 10) + ( (l' ( 1 0) "2. + 

iGU.) ) ' : II L- (• . • >) / a(2) +cf >* 2-2. • cs - ( 2 ) m.r ) / (T c + fUo) ) 

GP ( 3 ) = b.' I ; ) ■• ( . ( ; • ; ' 2 . ■ <; ( 1 3 ) ‘ ( VL--20 . ) - -.LOG I 0 ( P ( 1 1 ) +0 ( 2 ) /6 ( 3 ) ) ) + 

ic.i/ (.j( 15) . (• ( 14) nc n 



b‘f 

bb 

D6 

57 

5d 

59 

60 

6i 

6a 

65 

64 

65 

66 

67 

60 

69 

70 



17A 



Gin>(3) « oPP(a) -fGfM5) •i'^a-OPCa) ^x 2-a. + (B( 15) s-C2^ + cVD( 15) -i (GP(3) - 

lOrMa) )^ca)/(u( 14) 4 TO 

Cv ( 4 ) = OP ( ^ ) -^- ( j ( 17 ) > rB ( ,U ) *4 TK * 5 . ^ 0 ( 10 ) +T K + + 2 ) ^ C?0 G ( 2 ) /G ( 4 ) 

Gpr (U) = GPP (a) }-a,4GP(a) MCP(4)-GfM2) )H (2.+oaO) 4G,+0(i0)4TK)4C2 

1 > r i: *• 0 ( cO / 0 ( ) 

00 To 

b G.'ID (~«(2)+A(3) <.?,»T)-5/G() ) 

G,M ? ) = ( .\ ( b ) -2 . f t, ( 6 ) < 1 ^ 3 . - A ( 7 ) AT t *2 ) *S/G ( 2 ) 

Gr-!-'(2) = 1-2. f/\(6)+tic ^A(V) t-T)4b'!-+2/G(2) 

G;^13) - -0(2.) a(-A 19) +2c aA (iO ) aT- 3. tA ( 11 ) ATaa2-.'-4.»a ( 12) >Ta* 3) +S 
Gi’P ( b)= 2 . Alp (o) (5) A (2,* A ( J 0 )-6.*A( U ) AT K12. A A( 12) f T 1-A2 ) >Sa + 2 

GP ( 4 ) = ( A ( 1 4 ) “2 . t A ( I b ) + 3 . ( 1 ) #T-+.»2 ) *S/G ( 4 ) 

G,-P(4) n (-2.»Aab)-!-6. >A(lo)-?T)+StA2/G(4) 

6 Uu 7 1=1,4 
V G(i) = G( l)/GIIiP(I) 

PcTUR!) 

Kuo 
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ADAi-iS 

i 

. ^ 

a 

b 

6 

7 

8 
9 

10 

U 

lb 

14 

lb 

16 

17 

18 

19 

20 
21 
22 

23 

24 
20 
26 

27 

28 

29 

30 

31 

32 

33 

34 
3b 
3b 
31 

38 

39 
4 0 
41 
4 2 
4 3 
44 

46 

47 



SUl-ROU I I : JC. A DA ( ! . » H r X » I b.CI r Y » 0 Y r I ( IVE X > K ) 

D.;bt>LD PRECiSiO'i f’r Y; DYLLf DYLLLf L‘XL/DYiDYR,C2rC3rC4 rX^H 
Di.’-r.RSlOX' ^ n D) r OU lb) f P ( lb) . YR (i:.) ^D'fLLLI lb) rD'fLL( lb) rOYLUb) 

1 { lb) / C2{ ib) /C3(ib) f C4 Ub) 

1 IF( J5t:T*Gf ,C1) GO TO 6 

2 lr'(lij[.i^X«u.D*0) Go lO 4 

3 .< 2 

GO To b 

4 K r- i 

b C,,LL F(XpYfDY) 

b G-o To U0^7f8r9?ll) pK 
7 K = 3 
Go TO 10 
6 K = 4 
GO TO 10 
9 K r: b 

10 Dv 1001 lr.X9^ 

1001 P(.U = Y(J. ) »■ (H/2.0D0)-'Jc|jY( I ) 

Call F(x+:i/2*0D0rpfC2) 

Oj 1002 Ii^lWi 

1002 P(U ~ Yd) 4 (M/2.000)^C2(I) 

C/.LL F(X+h/2»COOrpfC3) 

Do 1003 Ii:ifN 

1003 P(i) = Yd) H* [-UC3(I) 

C.O.L F(X + hPp»C4) 

Do 1004 lrj,n 

1004 Y.iCl) ^ Yd) ^ (M/6«0O0)*(uYd)+2.QD0*:^C2d)*J-2,aD0:i=C3(I)4C4(I) ) 

Call F < X'hldUrDYR) 

Go 10 12 

11 Do llGl Irl/N 

1101 Pd) - Yd) > (n/24.0DO)^(bb.ODO-^<DY(I)-b9.0QO^OyL(I)-!-37,OQOiOVLLd 
l)-9.(rOO^'DrLLLd ) ) 

Call f(xhdpm;yr) 

d02 l-lpt'l 

1102 Y.dvl) ^ Yd) + (ll/24,ODO)-< :(9.0O0iOYRd)^19,OD0^nY(I)--b.0D0*DYLd) 

1 »• L.YLLd) ) 

C;.LL r(X«-|{; YRrUYR) 

12 X r. XiH 

Do UGl I'^XfU 

YU) - y;mi) 

D/LLi.d) ~ CYLL(I) 

DfLl.(l) = DYU(l) 
un.d) = OY(i) 

12 ji Old) “ oy:;(I) 

KATURi) 

r.iO 



1 

2 

:3 

4 

b 

6 

7 

O 

9 

10 

11 

12 

lb 

14 

15 

.16 

17 

16 

19 

2U 

21 

22 

25 

24 

2b 

26 

27 

28 

29 

50 

51 

52 

55 

54 

55 

56 

57 

56 

59 

40 

41 

42 

45 

44 

45 

4b 

47 

4 b 

49 

bO 

bl 

b2 

b5 
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S J r*v. 0 15 T I 0 ( » I ^ !n'*AX t H t Ai JGLc / Y < r~‘^ 5 Lf J t P[s Tf F ? f lEV » 5 r P ? 2 ) 

C b iPli. A ii f ,\LI- A y t Y » / OAi y G I GP ^ ? f 1 ( il r 0 V 1 1 )T / 2 V » P 

iy^PG0-»:jlroG/M5/DD 
L^GIC;.L .AjTPy.-.Gi.LtK 

CO, r 06 /i'L\ 4 /L-LPAiA/ 0 /PLK 5 /fi:;.WKT 01 / KMAK r KMAX 2 /GLK 7 / ALF AR ? ALP A I f 
|!v y O' '£GA/..'i.y.a/60TL (50 1 

uhAcubU) i £ (Oroy 5bl ) »L*V ( iriOvAX) fP (5 / 5/551) r KOiJ} i ( 35 1 ) / 0(14/701 ) / 
J..<( 5 ) ,F-(L 3 ) / r (py K/,A.< 2 ) / 0 .UPAR( 7 y 70 ,l) y YlNTvb) /OYiirr(G) 

£ .\ r L ‘ . I • ! A G 5 V f . o' c 0 
KPlOi - (K 10 T/ 2 )+l 
rjj t> I rb'(<=:i ? ITKMAX 
A.:-LPii\‘ ~ 90. 

c ESiAb... Ibr; uari AL COP'uITlOr.S A f liOGi'JDARY LAYPR EDGE 
IF (r.O «E( 14) }. v,Am: (orlOu) irEH/EVdTER) 
dlxPAP (/’ EV('»l2r») r;'/2Ll6,S) 



loO 

21 

22 

25 

24 

25 

26 
1 
2 

l04 



(21/22/25,24/25/26 ) t KLV 
= REAu(lV( ITER) ) 

“ AIiV,AG(LV'U rPR) ) 



REAL (F:V ( ITER) ) 
Aj..mA0(lV( ITER) ) 

RcAL(EV (ITER) ) 
'Au(lY(1TER) ) 



Go 1 0 
AiJ-AR = 

0 'LGA “ 

Oj To 1 
Ai.FAR = 

AuFAi = 

Go To 1 
AU'AR = 

R L ^ A 1 
GO To 1 
O..EGA = REAL(EV(I1ER) ) 

Re aIPAG(EV( ITEiO ) 

Go To 1 

OaiGA RE/aJEVdTER) ) 

Ad AI - AXMAG(EVdTLH) ) 

Cn) To 1 

R£ - .TEAL (r.Vd ter ) ) 

ALFAI = AI.nCdPA ( ITER) ) 

AuKA C..FLa{ALFAI;/ALP;U) 

IP (ALFAR ,LT. 0,0 aC5d O'.EGA ,LT, 0.0 -OK, f^E .LT. C.O) RETURN 9 
GA-T;-A ~ CSTThT ( ALF/\? + 2 M ,0 / 1 o 0 ) ( ALF 8LPAR ( 2 / KTOT ) --OPEGA ) ) 

b'GNA " CS jRvT( ALP A ^ ,C / I , ) f'RE U4U d ALFA njLPAR ( 2 aaoU ) -OPEGA ) ) 
If (L’OTEdL)) VKll'E ( 6 ; i iVi ) GAP.PA/SIoPA 
FoKPaT (» 0A\E,,A -» y2Fi:>,.Ty lOAdbJGf’A -*/2rl6,6/) 

IF (i<EAL(uA i:A) .LT, 0,0) fOA --GAP 1A 
IF (RdAu (bdPA) ,LT. 0,0) 5rC5\ - ‘-SIG.'d 

[ F ( ROl id 26 ) , /\i , 1 .'OT r 1 5<j ) ) GE» *" f 0 (10 c K I'O T ) •^ C ( / / KT01 ) ) ^ ( 

I Jlkar ( ^ / r To i ) '-C .•ega/.^- I r- a ) / c (i / l ror ) 

IF (,mOT„ ,;OTK(;:o) .c/l. rAPicdOO)) C5 r. (OrO/0,0) 



T.a ) 
roD 



2(L/i/Kr 
2 (A / I / Af 
/ ( 5 / j .♦ IN'* f oT ) 
Ed/p/Rprof) 
2 ( 2 y 2 y P lE; r ) 

2 ( G » • / »• r*’ I ) 

2U/C/''-Pr.d) 
2 ( 2 / / M I j I ) 
2(5/ 5.^P^■d ) 



= ( 1 « 0 / 0 , (• ) 

- ( 1 . 0 / u , u ) 

" ( 1) , u / 1 r U d 1 iF IA FiOi., 
= "*(i..O / I .0) :-Ai_FA 
~ , 0 d , 0 i i^GA.'G.’A 

“ mL/ .d K.l >jU/ (SIGF.A J 

-- i 0 , 0 y 0 , li ) 

- ( 0 dd »1 , ••; ) 

“ ( 1 , 0 / 0 , (• ) 



’ (EIGPA'i ^2-ALf A + ^2) 



;dAEFA*-^2 ) 
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56 

bV 

bo 

b9 

bO 

61 

6b 

64 

6b 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 
76 
79 
60 
81 
62 

03 

04 
8b 
86 
87 
06 
69 

90 

91 

92 

93 
04 
9b 
9b 
9/ 
90 
99 

100 
101 
102 
A 03 
lu4 
lub 
108 
ioV 



Do 3 \i~ 1 r i I 

/!( JM;,KPTO'f) = (6A.V.MA Jf l»KPTOT ) - (0 • » I . ) *2 ( J » 2 r KPTOT ) ) M ALPA + 

IGai A ) ( 0 . ; 1 . ) ( 4 . aLFA/o , ^ GAPFA+ ( ALFA+GA-^MA ) / ( /\Lr-v ^SI G;iA ) ) 

3; ^PTOT )/KE 

Z ( Jf fwKi^fvjT) = “-((Owl.) ^aLPA>CA’ ( Jr 1 f :<PT0T) + ( ALFA vGAFN’.A) + 

I2(wir 2# Kf^rOT ) ) iALr (ALFA f-oA’GGl ) ( J r 3 TOT ) / ( ALFAf SXGisA ) 

3 Z(JrorAPT^I) = -SlOFA-v:(J.-3::<; TOT) 

IF U;0lL ( iA ) ) UKlTi: (GrlOb) ( ( 7 ( J r I p Kf^TOT ) r I ^ ) f I r N ) 

li;b Fu?<F.Ar (' IrariAL values JLFOLl ooL. KiT£G’/( 12El0o4) ) 
ii - 1 

C IhTLO.<AT;: ijISTUKpAMCE LONS FOR AbSU'lEJ F, I GEL’VaI.UF p 0|;TH0N0KV,AL 12 XNG 
C iJlClSSAiF<" Oit VliEfcc specified 

Uu 4 r^Pp'^lrKProT 
Kp = r;PTO^^ I'-KPP 
IF (KP .LO, KHTOT) go TO 15 
Do 10 wL-irij 

K = 2 + Kp-i-l 
Do 9 I'*UM 

9 YlfJ(I) n Z(JrIrKP+l) 

CmUL uI FFEO ( STAhLO r DY IM r Y IN f r rW K ) 

Do lu 1=1 pM 
U' Z(JpI;KP) :: YINT(I) 

lb Call O.U!ICK(2p ITER rKPp ANGLEr II pPp KORTMp AGLMIN) 

4 C OK T I. SUE 

Ir (NOTE(i7)) LRITE (6pil7) AGLMM 

117 FokNY-T {* minimum angle detween SOUN vectors encountered in ORTHf 

IPPsOCESS =♦ f Fi0*6) 

Ip (AGLmIk .LT* 1.0) RETURN 9 
b IXEND - II - 1 

IF (f.OTEd?)) WRITE (6 p 107) IIEND 
lo7 FokMAT (' total no of ORTH REOD =»pI4) 

C evaluate SMI'SFFCTIOM OF B.C, AT TliE V;aLL 
3l = LL!bVR(4ri) 

D2 = (O.GpO.O) 

B3 = ( 1 . 0 p 0 • 0 ) 

X(i) = (IcOpO.O) 

DENO •. ^ (t;2tZ(3p6r l)+r53d(3p3pl) )J(Z(2p2d)+Bi'i=Z(2plpi))- 
l(j^:^^2(Ep6/i)+u3r/(2p3p3 ) ) 2 ( 3 i 2 # 1 ) +U 1 ^ 2 ( 3 p 1 p 1 ) ) 

X (2) = :< ( A ) t ( (E2'^2( 1 POP I ) >u3-!^Z(1 pop 1 ) ) ^ (2(3p2p 1 ) +D1 d(3f 1 p1))- 
l(o2+/(3p6pi)-J-D3-^2(3»3pj ) )^ ( 2 ( J. p 2 p 1 ) -^ R 1 ^ 2 ( 1 » 1 p 1 ) ) )/0FiJo;.) 

X ( j) = X ( i ) ^. ( (l2^2 (2 POP 1 ) +13.3 + 2 (2 p o p i ) ) + ( 2 ( i p 2 p 1) >Rl d. ( 1 p I p 1 ) ) - 
t X f I ) +R3 + 2 ( ipop 1 ) ) + (2(2< 2; I ) s-Blt-2 (2p 1 , 1 ) ) )/UENCM 

Dv lo 1 = 1 p •! 

16 Yud) = U ] )-5^2{ 1 p 1 p 1)4> (2) ^2 (2 p I p 1 ) +X(3) +2(3; I p 1) 

If (M‘v)ir ( Vv) ) v.oXTL (n;U2) ( Y ( I p 1 ) p I = 1 1 . 1 ) 

ii2 FvR’*- i (' r ( I p 1 ) 0* I i2-:io,‘0 

It (O o';( I ('• r 1) ) L.0ti'-2b) RETURN- 

r ( iTE-s) = I a / 1 ) /Y ( T d ) 

( i TEa) ) M E. J iter 

{ I «'. ■ E*. ( : i ) ) ./P X i ( 0 f X u O } i i ' » r ( d ■ ■( ) 



IF 

IF 

f,. 



, r 



6 C, 



rEO. 1) GO To 7 
• lTfr'RpCADS(f'( I ir.fO ) 



I 



r ‘dll.-. 
,AU‘E ‘ 

) p r. , : 

(T ? ur: 



) = ' P i. O • ■' / \ r * Of' C Do M ( 

ASoU ' E'j E T OC M Y a EL[_ 1 q US!.. 

> ‘ V p / ! I.Jv ♦ M.' V R ; R ) 



iD. f 



‘ ) 



d L 1 6 , E ) 
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IQh 
j..O 9 
110 
111 
I.U 

114 

11b 

lib 

U7 

1U3 

119 

I$i0 

.121 

122 

125 
124 
12b 

126 
127 

120 
129 

150 

151 

152 
155 

154 
15b 
136 
157 

150 
159 

140 

141 

142 
145 

144 

145 

146 
14 7 
140 
149 
IbO 

151 

152 

155 

154 

155 

156 



C 



IF (AUS( 1 ,.-REAL{cV ( IT.:K) ) /»^5AL(FV ( ITFRl'l) ) ) .GT* FPSUfj .CU, 

1A5S{ I ,”AI: AGvFv { I i'lP) )/Al.-:AG(eV { TT5KM. ) ) ) .GT. EPSLfs’) GO TO 6 
ir i nc\<) ) . 07 , cP-JLrn go ro 6 

ir {'!OiL{.:)) kETURN 
17 1-1 » 5 

l7 Xci) = X(1 )/CAl>5('i (471) ) 

Gj To 7 
o C0i,T IfTO'E 

IF (liOT..(2.l) ) ARIIL (6HCi) EV ( i U IF ) 7 f- ( N I TR ) 

r:ir5{^l) :: .FALSE. 

lei F.jfv^^.AT (/« ‘-AXli-U;. ALL.Ou'AjLE MO iTERATJONb FOR EIGENVALUE EXCf.EOElj 
1V» FEbr LSTH’.ATL oF E.V. = ' / 2E 1 4 . 3 ? ’ f ATTEMPT I LG TO SATISFY HOMO B 
2*C. WiT.j ,\b-1TTaNCE CFW2L14.b) 

IF ({iOr::{.?2)) uHIIl (6rU5) (EV( iter) iF( ITER) ? ITER“lf XTRMAX) 
ll5 roKF'A'f (15X r 'EV ( H t JO ’ 1 E.bXf » F ( I TER ) * / ( 4E16 * 8 ) ) 

GJ ro ( 51 52 f 55 7 54 7 55? 50) JviEV 
51 AlFAR “ REAL(EV{fil iR) ) 

O/'EGA = AlisAGCLVCFJTR) ) 

Kc.TU?^n 



52 AlFAR = RLAL(EV(fJI TR) ) 

ALIA I = AiFA6lEV(r:.lTR) ) 

RETURN 

55 AuFAR = UEAL(EV (MITIO ) 

RE = AlFiAG(EV(MITR) ) 

Return 

54 O...EGA = REAL(EV(MITR) ) 

RE AlMA(^(EV (faiK) ) 

RETURiN 

55 0,-iLGA = fiEAL(EV(NXTR) ) 

AuFAI = AlMAG(EV(faTR) ) 

return 

56 RE = HLALCuVCNITR) ) 

AcFAl = A.liiAG(£V(iJlR) ) 

RETURN 

E V aLUmT £ E 1 GEi J UNCT I0.\S 

/ IF (NOTEUL)) V.aITL (o?109) ( J ? X ( J ) i J--1 i N ) 

iu9 Format (» values of conjiiuic coeff at wall; 

1/(50Xf ’X( ’ I Ilf » ) iFElS.a) ) 

Call EO;,rC-j(X . Y? Zi IiF:;:P laXORTM) 

IF (./jlECiJ) ) v.KUl (o> 114) ( { 1 f XOR TH ( .U ) f 1 1 = 1 ' I ILND ) 

11^ ((» KCRTH( * 7 I5f ’ ) =»iI4)) 

M.aTF (6# 1 02) (KPf JiLL'AR { I f E:. KF-1} a Y ( I f kP ) i I = 1 / '. ) ? K.P = 1 # KI 



X( » ; Ilf * ) f2El5,8 



roT) 



l'j2 F jrj iAT (/Ijc 



r25X- 
I 1 1. X . 



• f ’KP' 

‘PRESSUIJP , 



»5X f NET A » f 6Xf » (RI.O+PN ( ) » f I’iXf »F» 
l4Xf »F« f JSXr ^TEFPE.RAfURr t/( I4f 



lloXf* Tl. pL.nAlUKE’ 

2F/f4flc.Elfj,5) ) 

W,U IF (Of U j) (hPrFLI-’AR ( : f L “1 ) I (C ''ESC Y ( I f XP) ) . I-l f ^’.) f KP::l »KPTCT) 
lib Fvl^.’LvT {» ...AI.B ihFiR COAFESPOr::: J.NG A ipL I TOOLS < / ( 1 4 f F 7 . ^ f 0E2 0 . 8 ) ) 
RElURiJ 



£m0 
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STAdEU 



1 

2 

b 

6 

7 



14 

lb 

16 

17 

18 

19 

20 
21 
22 
2b 
24 
23 
26 
27 
20 

29 

30 

31 

32 

33 

34 
3b 

36 

37 
30 



StJd'^O’jT I IL ST Ac E"Q ( Y ^ 0 Y ? K ) 

Coi- rLl-A i 7j>Y mVLFA^ t'lLtiSSC 
LjGlCAL u'OTF 

COf.'.roll /i^LK-Vbl.PA^ >C/i}LKb/rv>M^KT0T?KMAXrKflAX2/l5L.K7/ALFAK» ALFAI ' 
ir<G » OvlGA /,^LK6/LOTt. ( 30 ) 
in;,ZU:^l0 - uYC;) »YC-:) rbl.>-'A}<(7?701) ;C(l4r701) 

AuFa r: c .flx(alfa;naufai) 



a 


IF 


(;..OT;i 


( 


) CC'i'J “ 


1*0 




9 


If 


( .iJOT 


• 


..OT£(?e) ) 


coi. 


= 0.0 


10 


Ir 


( .Olf: 


( 


£o) ) DEN - 


1 .0 




11 


IF 


( ,;.OT 




iiOTC(E6) ) 


0 C_ i ; 


= 0.0 


12 


‘’*lE 


SSC r 




iiUF AK ( £ » h! ) • 




13 


ir 


(1X1 £ 


( 


jii)) CY(1): 


=•'(0. 


d » ) 4 



) M C ( 10 ? K ) “C ( 7 f K ) ) ^ ^ .vLESSC-Y ( 3 

1) ) 

IF (.EOT, note (30)) DY(i) - ( 0 • r 1 , ) +Y ( 2 ) 

Dm) - Y(5) 

If (FOTE(30)) DY(3) ^ Y(6)*C0N 
IF (.A-OT. i4OTE(30)) DY(3) = (0,0; 0*0) 

IF ({^OlF(dO)) DY('t) - (0 . ; 1 , ) t (“ALFA^^2 ^V;lLSSC-C ( 1 ;K) ( 1) + ( ALFA/R 

IE} ^ C ( 2 ^ A ) ^ ( " ( 0 * f 1 , ) * ( C ( o ; K ) “ ( ALFA ( 1 ^ X ) ) * ) -r Y ( 1 ) ^ ( C ( 4 ; K ) +C ( E ^ K ) 

2-? C ( 14 ; X ) ) . Y ( 2 ) -»Y ( b ) I ( ..LCSSC 4 DEN ( C ( G r K ) 4^Y ( 3 ) -5-C ( 7 ^ K ) <^Y ( 6 ) ) -s ( C ( 7 ? K ) 
y^DUl-C ( S / ) ) •«'FLPAF ( 3 > K ) ♦ Y ( 3 ) ) J'COM ) ) 

IF (,iJ0r. f;OTE(30)) DY(4) - ( 0 v r 1 . ) ( -AlFA + ^2'^LLES5C*C (1 ; K) t y ( 1 ) •}• 

1 ALFA/i-.E-':C 1 2 » K ) ♦ ( ( -Y ( 5 ) +Y ( 2 ) •+ ( 2 . ( 3 ? K ) +C ( 14 f iO /2 • ) ) /C ( 1 ; K ) -{' 

2(ut fl*U\i.FA + + 2*C(lF.<;)*Y(l) ) ) 

IF (iiOTE(oa)) DY(b) ~ ALFA*REt(BLPAR(3fK)*Y(l) ^(0,fl.)«(Y(2):^^ 
1 v/LESSC+C( 1 #K) + Y (4 ) ) ) ■^C( L r K ) /C ( 2 ; K ) v ( ALFA ( i ; K ) ) .'^•+2i<( (0 c ) 4C (4*K 

2 ) A Y ( J ) 4 Y ( a ) C ( 1 U r K ) -tWI-LSSC t Y ( 3 ) -^ DLi J-^CON ) +C ( 8 ? K ) *t Y ( b ) -f ( C ( 9 ; K ) 

34^^ ( 3 ) C ( L, I ,N ) 4= Y ( 6 ) 4 vil .PAF ( 3 ? lO ) ^ CON 

If C.iNOCe i;0lE(30)) [»Y(b) = ALFA Y ( 1 ) *bLPAR ( 3 / K ) 4 ( 0 , / I . ) * ( 

1 Y ( 2 ) ^ v.LMS::;^ i-C ( 1 ; K ) 4- Y C ( 1 f K ) 4 s ^/C ( 2 ? K ) + ALF A f ( C ( 1 f K ) ^'^2^^ Y ( 2 ) 

<14 (0 . ; \ . ) =4'Y ( I ) ^(C (^/ FO -C ( l'+rK)/2. ) i ^Y (b) :h (C(0f K)-C ( 14 ; K ) ) 4 Y ( 2 ) 4- 
3(C (3;K) /(■:. (4/K) -C (o»r.) )-C ( 14^L) ) 4C ( ;K) /2. 

If (NOTE(^O) ) LY (6) ( ALFA4^HL-*.M Y (1 )"^FLPAK (6?K)4 ( 0, f 1 . ) •4vvLESbC^Y( 3 

1) )*C(U;K) 4Y (3) M (Ml.FA^-^CdrK) ) 4=:^ 2-i-c ( 1 2 / N ) ) +C ( 13 r K ) :<Y ( 6 ) ) tcOil 
IF (.I.OT. hOTE( 30)) LY(6) = (O.OrO.O) 

RETURN 
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CCNFCil 



1 




r Iij;^ F.C'UrCii ( X i Y . Z f I 1 L':D> r- r KORTH) 


2 




I-,FLlClr COMPLEX (F-2) 






c ;■ . / a L k 5 / i ; i ■ ; > E r 0 T » K -l A X » K X ? 


4 




Dit r.;;EJO:i x(U) J Y (:.;»W-'.AX2) r Z ( N i M ' KMAX2 ) fQCJ) p KCRTlI ( KMAX2 ) . 


b 




lP(ljfi,.KVAx2) 


6 




li - 


7 




KplOT - (KT0T/2)H 


6 




D'j 1 .i" 1 p 1 J 


9 


1 


Old) ~ \ id) 


lu 




Od ') K“IpKPTOT 


J i 




Ir ( I I .1-E. 0) GO ro 3 


.12 




IF (i; pEE. KORUi(Ii)) GO TO 3 


13 




KTE.'-'P = KORTIKU) 


l-^ 




Ou 2- L=1 p;i 


lb 




C(L) - (0. Op 0,0) 


16 




Do 2 J=lnM 


17 


2 


0(1.) - 0(L) + P(l.p Jp|<Tr.MP) *X( J) 


10 




UJ 5 L=1p:i 


19 


b 


X(i.) " C. (L) 


20 




I J = 1 1 “ 1 


21 


5 


Do 4 I*" Ip 1*1 


22 




Y(1pK) = (O.OpO.O) 


23 




DO 4 J~1pi4 


24 


4 


Ydp.O = Y(IpK) + Z(JpIpK)+X(J) 


2b 




UElORiJ 


26 




E,pO 
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osorth 



1 




SJI.ROWT r;,. CbOMH(Z»PfK) 






CO: I'LLX I-' j > T iZ 






L.j'jIcAL. .'JUTE 






C •? Cl ; /RLKb/;-i ' ' !• TOT » Kf-'. AX > K.MAX2/bLK 3/M ATE ( 30 ) 


5 




Oi.- Ei.SiO.i S(3-3) »V, (3i3) . T(3>6) » Z ( iM , f'. # KMAX2 ) f P UJ f N , KMAX2 ) 


o 


1 


OC lo 0=1 »M 


7 




Oj a i = l>,l 


a 


2 


T(Jfl) = (O.OfO.C) 


9 




IF (M .EO. 1) 00 TO 6 


10 




j:.top = j ~ 1 


11 


3 


O.J b L=1»J'3T0P 


12 




S ( « L ) = ( 0 . 0 . 0 , 0 ) 


13 




Dj 4 1 = 1 M'l 


1^^ 


4 


Su»i.) = S(J>U + CCMJG(2(L» I rK) ) -/( Jf I »K) 


13 




Uv/ I“1 f 1*1 


1.6 


b 


T(J»I) = T(J»I) “ S( JfL)*Z(E» IfK) 


17 


6 


00 7 1=1 ;M 


16 


7 


T ( u » 1 ) = T ( J » I ) ;: ( J » 1 f K ) 


19 




vn J » J ) = 0.0 


20 




DO 6 1 = 1,;., 


21 


8 


= -.-.(JfJ) + COI!JG(T(J,I))^T(J,I) 


22 




w(j,j) = sonT(w(j,j)) 


23 




IF (.•,(-•, J) .LE. 1.0E-2Q) RETURN 






DO y I-l»:-i 


25 


9 


ZUI'lfK) = T(J>I)/.<(Jf J) 


26 




IF (note (2)) GO TO 16 


27 




DO .15 L=i,n 


26 


lO 


IF (1.— -') U,13,J.4 


29 


11 


P(L»J,R) = (0.0, 0,0) 


30 




DO 12 LUL,JST0P 


31 


12 


P(L»J,K) = P(L,J,K)-S(J,LI)<=P(L,L(,K)/W(J»J) 


32 




GO TO IS 


30 


i3 


p(L,j,K) = c;:plx(i.o/w(j» J) ,0.0) 


3^ 




G J To IS 


35 


14 


P(L»J,X) = (0.0, 0.0) 


36 


lb 


Cj.'.Tii.UE 


37 


16 


C.;j;TII;Ue 


38 




return 


39 




EnU 
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tVLST 

1 

2 

,5 

4 

5 
b 
7 
b 
9 

10 

11 

12 

15 

14 

I'b 

lb 

17 

lb 

19 

20 
21 
22 
25 

24 

25 

26 
27 
2b 
29 

50 

51 

52 

55 
54 
5o 

56 

57 
56 

59 
4 J 

41 

42 
4 5 
44 
4b 
4b 
47 
46 
49 

60 
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s I f r. v-:s t ( p ^ i ti-ik ax ^ vj f i iua\ , :.I T2 ? v.ull^u ) 

Oi;:LliSiC‘'i K( 1 TKMAX) iTH.-.iAX) 

I :l LICIT CO“.f'’LLX (A“‘H/f’) 

Ln^GI CAL .’'iOl I" t : 'L‘LI_l.K ►' avg 
/ uL:\ /bOTt:(50) 
il If- UTLT . 5) CO rO 6 

A;b :: .FALSL, 
r-i^^LLiCR - t FALSE# 

IJiTK = 1 
Oj 6 lTRn:l,5 

Ir (CA'6S(F ( TTiO ) .LEc C AQS ( F ( L i TfO ) ) EITR ~ XTR 

6 COlTIhUE 
Gj To / 

b IF (CA6S(F( nAA-:;) ) .LE. C ASS ( F ( lA I TF ) ) ) ;llTf^ = ITeR 

IK (CAES(Fa Flk) ) ,LE, b*uL-5) AVG = .FaUSE, 

ir (ITER ,iiL, (ar; ) avg = .Tf.UE# 

If ( i :0T E ( 2 5 ) ) L H IT E ( 6 > 1 (; 4 ) li I i la C A : J S ( F ( r,' I T R ) ) 

i04 FOivl'lAT (* CAr:.S(F( » » 12? » ) ) ::»?E16.8) 

7 2i = Air;Ao(F(ilER-^2)-^C0!jJG(F(lTER-i))) 

r Ali:AS(F( lT£K)^COilJG(F( iTER~2) ) ) 

2j - Al.-AS(F( lTER-i)+COLJG(F(XTER) ) ) 

S n Ab5(RE.,L(LV(ITER) )/nI AG ( EV ( I TER ) ) ) 

If (S *Gr, l.GEo .OR. S .LT. i.OE-’S) GO TO 9 

If' ( (2l + -^2'i-22<^=f2«*Z5'^*2) .LE. I.E-'IS .Arn«>*UOT. N.UlIJER) GO TO M 
9 E,/(irEK + l) = (ZltEv(ITErJ+Z24.EV(ITEF?~l) fZ5<-EV(ITER-2) )/(21-^Z2}-Z5) 
If (AVG #OR. nEU.i£0.5) EvdTERM) = ( E V ( ITER+ 1 ) ^EV ( N UR ) ) 'M # 6 f . 0 ) 
MJLLER = .FALSE, 

RuUiRU 

a IF (HER #E0. 5) GO TO 2 

1 IK (CAbS(F(ZTEr<)/FnTER-I) ) .LE. 10.0) GO TO 2 
EUITERU) ~ EV(ITER) + (CV(ITEP)-EV(lTFR-l) )/2#0 
IK (;iOlE(27)) AKITu (6?. 105) 

lu5 FjRFaI (» bTFFEREiiCE oETV.EFJJ F(ITER)S ToO LARGE To USE F.UlLER,») 

Rf HURM 

2 lU - EVdlER) - EVdTER-i) 
riiKl = EVdTER-d) EV(lTnR-2) 

El - hi/HlMl 

U1 ~ 1*0 I- El 

Gr F ( I rdR-2) =^EU^2 - F ( I TER-1 ) ^ C I L F ( I TER ) »:‘ ( f I + fj I ) 

RAU r CSORT (GUI 2--4. ^F( I TEK) ^DUE. I ^MFd TER-2 ) E 1 -F (1 TER- 1 ) =«^D I ; F (I T 
ILR))) 

u.^iiOKi ^01 + KAO 
G 1 . f A U - o 1 "" F' /» C) 

Ik (CAbS(0.:L'0:i) .LT. C ALS ( G I : :R AD ) ) CENO I = GlFRAD 
J:- (C/^6S(l:.LK G.K) #Ll-. l.CL^‘20) GO TO 4 

EJ>1 := -2.0^K (UEK) +di/uElo;; 

5 Mil ), r. FUn^MI 

n . (1 ( t -d.l ) = LV ( I FcR) V mPl 
J- (aJ’KC:/)) v.MTE (CdOl) EVdTLf<-^l) 
lL.i F..RUU (» ..LuLLi'.S ...CT‘K(; IM.OUIF^EO FOR f.iIS IlERATiOM. EVdlEKal) n 

1’ 16. a) 

;.:^LLfK« = . n:'Jh . 

Rl ) (»f;k 
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5^^ 


'» Eii'l r (1 


bb 


IF 


bo 


l^j^ F>.Ki ;;,T ( • 


57 


Gu ( 0 ^ 


bO 


f-„0 



. 0 ' 0 . L' ) 

;y)) v.KiT(£ 

i-IPJ bET 



(6f 

touAL ro i . 0 • ) 
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ORliiCK 



1 




SuiiKOUTl;:._ Oirri.CIUZ* nEK»Kr ANOLEr IJ / P t KORTlI / AGLmIi 


. ^ 




Cj.-VLcX P,Z,S 


,3 




L.i;0 1 C aE f j J f E 






Cg:-,; ' ofi /o'..;vS/.i . : ; » PTOT < k :aX > kmax2/elk0/mote ( 30 ) 


b 




0 it iE' .S 1 0: 1 / ( IJ , I , :C aa2 ) , P ( U h'J . K.'-iAX^: ) f KGR HI ( KMAX2 ) ? : 


6 




if (k «E0. 1) (rU TO 6 


7 




il- (AiA-I.E .0!.:, yo.u) GO TO 6 


ti 


C ilf'E BEli.EL'Il SOl.i; VECTORS 


9 




!)0 4 ,j—J. ?H 


10 




SlJrJ) :: 10.0,0.0) 


Xi 




Ou 3 1 = 1 . 1 


12 


3 


S(JfJ) = S(JrJ) + COMJG(Z(J»l fK) )+E(J»I»K) 


ib 


4 


S(vjfj) r (S( J>J) ) 


14 




OvJ 9 L=2/n 


lb 




JO = L--1 


16 




OJ j 9 0“*1 , 00 


17 




S(J.L) = (u,0rU,O) 


13 




00 S I 1 f Ai 


19 


b 


S(J»E) = S(JrL) + C0!10G{Z(J»Ifl<) ) >2(L..I,K) 


20 




G.\W'iA = A'.'jS(CABS(S( JrL)/(S(J. J)+5(L.-L) ) ) ) 


21 




GaI-O-.A = GA'.A+Sy. 20378 


22 




AGLfiltJ = AnIIU (AGLNON- GAMMA) 


23 




!r ((..'0TE(2'') .ago, I TEK .EO. 1) V.RITE (S'lOl) J»L 


24 


lul 


FOIO- Al {' AMGLE bETV.EEli SOLM VECTOR* rI2»' AMD’ » 12 


25 




1* 1S'»F<.i.4,' LLOREES.') 


26 




Ip (gamma «LE. argue ) CO TO 6 


27 


9 


COMTIijUE 


23 




KETURR 


29 


6 


C„LL G30RTH(E.P.K) 


30 




K^jKTh(U) = K 


31 




11 - I 1+1 


32 




RETURN 


33 




END 
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sjL)f:ojri .‘c QiFRjvd- h:)Y^ Y Mi^K) 

C.j.:!’LFa 'JYi Y f f’ rC2 f C3f C4 

Co *•: .Ci,' / JLK5/,\ » i'i f K 1 OT ^ KV;AX ^ KMAX2 

0 i ’•.L.f ] 0 I 0 Y ( M ) /V ( i . ) f C 2 ( b ) f C 3 ( b ) f C 4 ( 6 ) f P ( 6 ) 

C,,LL r(Y,,jYFK) 

Do I 1-1/M 

1 P{\) ^ YdJ + (H/2.0)*DYa) 

C;,LL r(PiC2F<-l) 

Dj 2 }=1 fM 

2 P(l> Y( i) + (M/2.Q) iC2(Y ) 

Call F(pFCjFK-i) 

D .) 3 1=1? : 

3 fMI) = Yd) + MtCo(I) 

C^LL K(PfC4fK-L) 

Oo 4 1 = 1 F Pi 

4 YU) = r(X) + (H/6»0)ic(oY(I)+2,^:^C2( l)+2,+c3(l)-fC4a) ) 
He 1 L'pjvl 

L.iO 
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The Coi fipulre r Pro gr ain in Opera tion - 
Selection of Nun ieri c al Param et ers 

The numerical program developed to solve the stability pro- 
blem outlined in Chapters II and III is written in FORTRAII V lan- 
guage and evaluated on a Univac 1103 digital computer* Just as 
discussion of the complete problem has been divided into tiw 
categories (mean and disturbance) j so too is the computer program: 

(1) Solution of the mean flow equations in real, 
double-precision arithmetic, and then cal- 
culation of the. variable coefficients required 
for solving the disturbance equations; 

(2) Solution of the eigenvalue problem in. complex, 
single-precision arithmetic and subsequent 
manipulation of the eigenfunctions to compute 
Reynolds stresses, dissipation, etc. 

Operation of all or part of these two elements is regulated 
by a series of input logical variables. 

In addition to the parameters (or eigenvalues) vahich expli- 
citly characterize the forn’iulated stability problem (a^, a^. , o) , 

L\ 1 ^ 

Re, and of course and M) , there are also adjustable para- 

meters which implicitly characterize its numerical solution (e.g., 
step size, orthonormalization erriterion, specification of boundary 
layer edge, etc.). Since the latter govern the ability to 
accurately and efficient.ly determine the former, discussion of the 
program's operation \;ill be developed around the optimum selection 
of these numerical parametcr.s. 
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To better illustrate the difficulties associated with 

making such a selection, numerical examples are calculated for each 

parameter for the saxae typical "eigenvalue^* with fixed: 

T - 90®F, T := 60°F, H = 0. = 0, and Re = 7000. 

V7 * CO ' I ’ 

Table E.l demonstrates the sort of results obtained for the same 
fixed parameters and v;ith: h, “ .0175, Q ~ 88'', p =10, and 

n. = 5.7925 (u^ = .999020). 





L) XIO^ 

CO 


CoTUTTients 


.123764 


3.72957 


Including fluctuating temperature 
field but without density 
fluctuations . 


.121905 


3.66240 


Without fluctuating temperature 
field but v;ith all irican fluid 
properties variable in dis- 
turbance equations. 


.122347 


3.68283 


Including all mean and fluctuating 
quantities, equations (2.37 - 2.^0) 


.118809 


3.56218 


Without fluctuating temperature 
field but with variable mean 
viscosity only (VJaz.?:an, et.al 
equation (1.2) [3 - 6]). 


.119928 


3.60135 


V/azzanVs equation using property 
variation described by Kaups and 
Sraith [7] (see Figure B.7). 



TABLE E.l The Influence of Included Disturbance Quantities and 

Chosen Property Variation on an Eigenvalue Calculation. 
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E.l Boundary Layer Edge, n and Hr 
max 6 

T\7o distinct length scales must be specified to accurately 
solve the eigenvalue problem: n , the finite value of the 
independent variable, p, at which the mean flow's asymptotic 
boundai*y conditions (eqns* (2,21a)) are assumed to be approximately 
satisfied; and » arbitrarily selected boundary layer 

edge V7hcre asymptotic disturbance flow conditions are satisfied. 

For computational purposes, there is no reason to suggest that the 
two need be the same* Wiile the former is associated \;ith obtain- 
ing accurate mean velocity and temperature profiles, the latter 
determines \;here these profiles may be truncated for the eigenvalue 
problem without appreciably degrading the accuracy of its solution. 

The practical capability for determining the proper 

however, is contingent upon utili^^ation of an efficient raeans for 

calculating the initially unspecified, mean-flov;, wall boundary 

conditions, u'(0) - F"(0) and H.’(O). The effectiveness of the 

particular laethod used in this investigation [18] for finding these 

conditions is demonstrated by the example in Table E.2 vjhere 

T - 90^F, T = bO'^F, H = 0, h = .02, p = 10. 
w > CO ’ m max 






C,r 











'~^ p- 
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Iteration n F“(0) H'(0) 



Guess 




1.0 


-1.0 


1 


2 


0.63890603 


-0.80589459 


2 


10 


0.45716795 


-0.71569664 


3 


10 


0.46637891 


-0.71328643 


4 


10 


0.46639972 


-0.71321361 


5 


10 


0.46639971 


-0.71321361 


6 


10 


0.46639971 


-0.71321361 



TABLE E,2 Convergence History of Unknounn 2-iean Flow Boundar>^ 
Conditions, F’’(0) and H'(0)* 

With 0*55 seconds required for the first iteration, approximately 
2,78 seconds/iteration thereafter, and six iterations, both con- 
vergence to values of F‘‘(0) and H*(0) with better than eight-place 
accuracy (for the specified fluid property variation, of course), and 
calculation and allocation to storage of all tnean-flow coefficients 
needed for solution of the disturbance equation can be accomplished 
in less than twenty seconds (see Table E, 3), 

Final selection of p is determined by observing the value 
for which an increase in this parameter fails to make a discernable 
change in F’*(0) and H’(0). It is reasoned that values of n 

ID 3.^ 

greater than this one will not alter the mean flow profiles (thus, 
the variable coefficients for the disturbance equations) and so 
must represent an unnecessary expenditure of conputational effort 



and tine. 
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The results of such an n search for the mean. V7all 

max 

temperature and velocity gradients considered in this stud^^' are 
presented in Table E.3. In all cases, F”(0) and H*(0) v;ere un-- 
changed for n >10. Thus, noting that the time required to 
solve the mean flow equations a.nd compute disturbance equation 
coefficients varies only slightly over the n range considered, 
accuracy rather than economy governed the selection: 



n 

max 



10 . 



This value ensures satisfaction of asymptotic conditions (eqn. 

8 

(3.4)) to order (10 ‘ ). Due to its ivaportance in calculating the 

eigenvalue, the influence of p on is also displayed in 

Table E.3. As a basis of comparison, it may be noted that the 

value of F”(0) for the unheated, flat plate case agrees well with 

Howarth*s value of F”(0) = .33206 [9, p. 129]. 

The primary considerations for selecting p^ were founded in 

attainment of an accurate eigensolution and reduction of running 

time and storage requirements (known a priori to be relatively 

large v/hen eigenfunction recovery is required [20]). Although the 

time restriction could be relaxed somewhat if the inforraation x;ere 

needed badly enough, storage limitations (150,000 for the Univac 

8 

1108 computer) could not. To evaluate the effect on the eigenvalues 
of usiTig a ‘'truncated” mean flow solution, the example shown in 



Table F.4 was calculated with: 
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n 

max 


F"(0)-u'(0) 


H'(0) 


E 


(defined 
by Eqn. 


Time 

Required 












(sec) 


T = 


60®F 










6 


.33255559 




1.7152416 


.52436x10"^ 


13.8464 


7 


.33209439 




1.7203045 


.43557x10”^ 


15.3772 


8 


.33205907 




1.7207617 


.13211x10“^ 


16.7630 


9 


.33205738 




1.7207868 


.14727x10"^^ 


17.7298 


10 


.33205733 




1.7207876 


.60392x10“^^ 


18.8970 


11 


.33205733 




1.7207876 


-20 

.91101x10 


20.9158 


T = 

M 


90°F 










6 


.46680267 


.71341605 


1.5112973 


.19198x10“^ 


14.6256 


7 


.46642702 - 


.71322750 


1.5142917 


.12898x10"^ 


16.1574 


8 


.46640087 


.71321419 


1.5145365 


.31727x10"^° 


15.2932 


9 


.46639974 - 


.71321362 


1.5145486 


.28695x10"^^ 


17.0454 


10 


.46639971 


.71321361 


1.5145491 


.95693x10"^^ 


18.9308 


11 


.46639971 - 


.71321360 


1.5145491 


. 25249x10“^^ 


19.7002 


T 

w 


= 150''F 










6 


.74345286 - 


.78390563 


1.2098909 


.35870x10"^ 


14.5934 


7 


.74320905 


.78382177 


1.2111477 


.17598x10'® 


14.5190 


8 


.74319434 - 


.78381662 


1.2112361 


.31693x10'^^ 


15.8284 


9 


.74319379 


.78381642 


1.2112399 


.21000x10"^^ 


16.9238 


10 


.74319378 - 


.78581642 


1.2112400 


.68194x10*^® 


17.7392 


11 


.74319378 


.78381642 


1.2112400 


1 R 

.17341x10' 


19. 3116 


TABLE 


E.3 Influence 

and rir*. 

6 


of n on 

max 


computed values of F”(U) 


, H’(0), 
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T = 
w 


200"F 












. 6 


.96336141 


-.83923358 


1.0362211 


.11659x10' 


■6 


12.3970 


7 


.96318694 


-.83918394 


1.0369269 


.47597x10" 


■9 


13.5738 


8 


.96317728 


-.83918114 


1.0369723 


.71388x10' 


■12 


14 . 8100 


9 


.96317695 


-.83918104 


1.0369741 


.39403x10' 


•15 


16.5640 


10 


.96317694 


-.‘83918104 


1.0369741 


.68671x10' 


-18 


18.4704 


11 


.96317694 


-.83918104 


1.0369741 


.61905x10' 


-18 


18.5032 



T - 60"F, M = 0, h = .02 
® ’ ’ m 

TABLE E.3 Contlm:ed 



= G0°F, H = 0, = .02, Q = 45", Re^* == 7000, a^. = 0 (lower 

branch v.’ithout density fluctuations) . 



"6 


^6 


i? Steps 


// Ortho- 
normali- 
zations 


Time/ 

Ite.ration 


“r 


0) xlO^ 

CO 


7.00 


.999959 


351 


102 


11.1514 


.12"3762 


3.72956 


6.70 


.999903 


336 


102 


10.7356 


.123765 


3.72957 


5.80 


.999033 


291 


99 


9.5320 


.123817 


3.72989 


4.68 


. 990055 


235 


96 


7.6570 


.124466 


3.73420 


3.14 


.900150 


158 


91 


5.1472 


.134641 


3.83808 



TABLE E.4 Influence of the Arbitrarily Specified Boundary Layer 
Edge, ri^, on the eigenvalue computation. 

Note that consideration of additional, segments of the 
boundary layer does not appreciably change the number of ortho- 
nonnalizations required, suggesting that most of tiicm are required 
in the region near the wall. It does, however, require more 
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integration steps and so more time per iteration. Unfortunately, 
due to storage, and time lirrdtations , it is more difficult to as- 
cribe a definite value to (also the case for the disturbance 

0 

equation step size, h,) than for n , F.ven using the maxiinun 
allocated storage (h^ = ,02 and = 7.0), variations in the 
eigenvalue arc. still slightly greater than their indicated six 
place accuracy. 

Since even four place accuracy is adequate to clearly 
delineate between calculations with and without temperature 
fluctuations (see example in section 3,4), it is sufficient to 
specify 

= .999 

6 

as the position at v;hich the boundary layer may be truncated for 
stability calculations . 

Greater accuracy than that indicated in Table E.4 is 
restricted by both storage and the interdependence of p^ and h^^ in 
the prograra; v;hen h^ )nust be reduced, so too must p^ according to 
the relationsliip > p^/350. 
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E • 2 Step Si ze, h 

Perhaps no other parameter more critically governs the 
accuracy and efficiency vjith which a solution is numerically cal- 
culated than does the step size. If it were chosen too large, the 
fourth order Runge-Kutta and Adaras-HouJ ton integration schemes 

5 

(with truncation errors proportional to h ) would be basically 
inaccurate, rapid changes of the dependent variables (such as in the 
thermal boundary layer for large Pr) would be obscured, and (most 
importantly for this analysis) between consecutive steps, for a 
given aRe, linear independence of the solution vectors may be 
destroyed beyond accurate restoration by the or thonormalization 
procedure, independent of the specified (see section 3.4.3). If 
chosen too small, not only v/ould integration time and storage be 
prohibitively large, but error propagation, enhanced by the in- 
creased number of mathematical operations required, would contami- 
nate the solution sooner and thus require more frequent ortho- 
normalization \7ithin a given interval (see section 3.3), 

Thus, to solve the disturbance equations, the selected h 
must be very closely related to the specified values of aile and fi. 
Since accurate convergence is ensured only by giving the ortho- 
normalization procedure a chance to work effectively, as aRe 
increases, h should be decreased and il increased. In mapping out 
a stability contour, the point at udiich convergence difficulty can 
be anticipated is preceded by a limited range of smaller aRe in 
vjhich the number of Iterations required for convergence is greatly 
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increased and iiiore importantly, the rainimura angle encountered be- 
tween solution vectors drastically drops to less than 1^. At or 
be3'Cnd this poixit , convergence could be obtained only be decreasing 
the step size. It is interesting to note that similar difficulties 
v-rere encountered by Hack [24, p.2.S3] (in a scheme not using 
orthonomalization) , but at much smaller values of aRe, 

In contrast, for the numerical values calculated using the 
fourth order system of Wazzan, et,al [3], no sudi difficulty v;as 
found. One possible explanation can be advanced by considering 
the examples given at the beginning of section 3,4. For the sixth 
order system (that is, the problem as formulated in this study but 
without density fluctuations), the minimum angle encountered during 
the boundary layer integration was 36,40^, whereas for the fourth 
order system (Wazzan^s equations) it was 78.21^. Since the region 
in v;hich significant thermal fluctuations exist is much smaller 
than that of the velocity fluctuations, it can be expected that the 
step size accordingly must be smaller for the sixth order system 
to maintain the same degree of orthogonalit}^ as for the fourth. 

Since the Runge-Kutta integration scheme requires infor- 
mation at half-step intervals, the moan equations are integrated 
V7ith a step size, h , half that of the one used to solve the 
disturbance equations, h^, (i.e., Tliis permits the 

program to calc\ilate the required variable, mean-f low-dependent 
coefficients for the latter exactly for all points. 
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An example illustrating some of the points discussed above 
is given in Tables E.5 and E.6 for: = 90 °F, = 60“F, M 0, 
n == 10, 11. = .999, Q = Re.* = 7000, a^. =*= 0 ( lower branch 
v;ithout density fluctuations). 



h 

m 


u'(0) 


ll'(O) 




.01 


.46639971 


-.71321361 


1.5145491 


.02 


.46639971 


-.71321361 


1.5145491 


.03 


.46639972 


-.71321358 


1.5145490 


.04 


.46639973 


-.71321352 


1.5145489 


.05 


.46639977 


-.71321341 


1.5145489 


.05 


.46639982 


-.71321321 


1.5145487 



TABLE E.5 Influence of Step Size, on the Unspecified Wall 
Temperature and Velocity Gradients, u’(0) and K (0) . 



"d 


//Steps 


Time/ 

Iteration(sec) 


//Orthonor- 

malizations 


“r 


w xlO^ 

CO 


.01675 


347 


11.1042 


103 


.123819 


3.72996 


.02 


291 


9.5320 


99 


.123817 


3.72989 


.03 


194 


6.4768 


94 


.123812 


3.72967 


.04 


146 


5.0016 


78 


.123757 


3.72917 


.05 


117 


4.0280 


67 


.123774 


3.72834 


.06 


98 


3.3430 


59 


.123741 


3.72718 


TABLE E. 


6 Influence 


of ti*)e Step Size 


, on the 


Cennputed 





Eigenvalue . 
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From these calculations (aRe = 866.6), it appears that at least 
four place eigenvalue accuracy is obtainable by using < .03. 

Note that this choice also ensures eight place accuracy of the 
mean flow solution (for the specified property variation) . For 
results given in Chapter IV, a step size of - .02 or smaller was 
used. 

^ ^ Orthonorm a lizatio n Angle, C l, 

As previously indicated (section 3.3), the orthonormalization 
criterion used for this study is a specified uiininura angle per- 
mitted between complex solution vectors. (A comparison of relative 
magnitudes of the vectors also could have been used [39]). If 
this angle Vv^ere chosen too large (e.g., 90^ as for Wazzan, et.al. 
[3]), aside from an increased integration time, the extensive 
matrix manipulation required for the orthonormalization procedure 
could be expected to introduce errors in the solution [39 ] , 
Alternatively, if it were chosen too small, then linear independence 
of solution vectors could be so degraded that insufficient resolu- 
tion would remain to accurately reestablish an orthonormal basis. 

Since the parametric group aRe appears exponentially in 
the exact solution to the disturbance equations outside the boundai*}^ 
layer (see equations (3.11 - 3.15)) or inversely as a coefficient 
multiplying the higest order derivative term (see equations 
(1.2 - 1.4)), it can be anticipated that as aRe increases, so too 
should the rate at v.diich the solution vectors linear independence 
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degenerates. It folIov7S, then, first that for the specified vQ, 
the number of orthonormalizations required laust also increase. 
Secondly, as Q is increasing to accommodate the larger value of 
aRe, so too should the percentage of orthonormalized integration 
steps, ranging from zero v/hen Q = 0® to one hundred when 9. - 90® r. 
These observations are substantiated by Figure Fv.l and Table E.7. 

Care must be exercised to ensure that Q is chosen suffi- 
ciently largo that the resolution of solution vectors is not lost 
betv7een successive integration steps just before the orthonor- 
inalization criterion is violated. That is, if the smallest angle 
betv/een solution vectors at one step is greater than but less 
than 9 at the next, this ninimum angle must still be large enough 
to enable reformation of the orthonormal basis. RTien degeneration 
of orthogonality is sufficiently?’ rapid to require orthononaaliza- 
tion at consecutive steps, it can be expected that the minimum 
angle between solution vectors (and thus convergence capability) 
v/ould cease to be a function of 9 and become one of the step size 
instead, as indicated in Figures E.l and E.2. 






1 



Mininutn angle encountered in last eigenvalue iteration (degrees) 
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OrthononaallzaLion A.ir,]c,Q (dcKrcps) 



Fif'.urc E.l Influence of the Ortlionorral iz/ition Anej e, on tlic 
nvirnber of ortlionorinall irntions required ]>ct: iteraCiun and on the 
lalninun nrifjc encountered in the laat Ijoundary layer Intep.ration, 



percentage of steps orthonorrnali 



inimum Anglo Encountered in the last integration (degrees) 
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9 . 0 ^ 



8.0 



7.0 



6.0 



5.0 



4.0 



3.0 



2.0 




1 j j. 



0 



.02 



.04 



.06 



h 



d 



Fip.ufe F.2 Influence of the 5tc-p size on the minimum onplc 

encountered bf*lv;eoi' solution vectors on ihe 1/ist eipcnv.sluc 
iteration for n specified 0. 



Q 


“r 


(j xio'^ 

CO 


Number of 
OrthoTiorraaliza- 
tions required 


Minimum 

Angle 


5 


.179330 


3.70852 


44 


1.3547 


10 


.179330 


3.70S52 


50 


1.5390 


20 


.179330 


3.70852 


74 


1,6808 


30 


.179330 


3.70S52 


79 


1.6808 


40 


.179330 


3.70852 


83 


1.6808 


59 


.179330 


3.70S52 


88 


1,6808 


60 


.179330 


3.70S52 


93 


1,6808 


70 


.179330 


3.70852 


100 


1.7733 


80 


.179330 


3.70852 


112 


11.3345 


88 


.179330 


3.70853 


140 


12.3368 


90 


.179330 


3.70853 


146 


not computed 


where 


Re = lO'^', 




.04, - .999 (upper 


branch, r 0) 



TABLE E,7 Influence of Orthonormalxzatiou AnglCj $7, on Couputed 
Eigenvalue, 

On the basis of these calculations, it i^ould seem that al- 
though the ninimum angle between solution vectors is profoundly 
affected by the specified orthonormalization angle in some 
regions, the computed eigenvalue is not. Further, in the 
intermediate range of there is a region in which the number 
of required orthonormalizations is significantly less than for 
Q 90® and v;here this number varies only slightly v;ith 
Naturally, such a reduction translates as a significant savings 
in required program running tire. For example, for the Re = 7000 
(upper branch) calculation shoi;n in Figure E.I, tlie time required 
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to complete eigenvalue convergence (in seven iterations) is re- 
duced fron 36.9 to 34.4 seconds by specifying Q equal 50° rather 
than 88° j still V7ith almost six-placc accuracy. Thus, selection 
of Q apparently can be based primarily on computational economy; 
that is, it is chosen to reduce the number of computer manipulations 
and program storage requirements. 

Since these results are typical of others found during the 
course of this study, it is felt that the selection 



Q = 45° 

represents a good balance betvjeen accuracy and efficiency. To note 
in passing, this choice is also compatible with the results of 
Gersting and Jankov/sl:! [20], 

E . 4 Eigenvalue Est imates, a, to, and Re 

As indicated in section 3,3, three eigenvalue estimates are 
required to start either of the root-finding schemes used in this 
study (i.e., that of Wazzan, et.al, or Muller). To trace out a 
contour in parameter space, the first eigenvalue estimate, 
taken from the lo^/er branch of the stability curves of Wazzan, et.al, 
[4] near the rainimum critical Reynolds number (see for example 
Figure 4,4 ). This procedure is employed for two reasons. 

First, since the present investigation represents a departure from 
the Hazzan results due to the inclusion of a f 3 actuating temperature 
field and consideration of all fluid property variation, it \;ould 
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seem that these should give the best available initial estimate. 
Second, as the difficulty of solution increases v;ith increasing 
aRe (due to problems associated with degeneracy of the orthogonality 
of solution vectors, etc.)> eigenvalue on the new contour is 
chosen for which solution should be easiest (see also discussion by 
Mack [24, p. 280]). 

The second and third eigenvalue estimates are arbitrarily 

chosen as ev “ 1.01 ev and ev^ = 1,005 d- 1.02 ev , With 

Z X j IK X i 

the three initial eigenvalue estimates thus specified, the three 
corresponding boundary layer admittances, Y^, can be calculated. 



conditionally followed by that of Muller. That is to say, only 
v/hen all the following criteria are violated will the latter be 
utilized: 



The plane fitting technique of V7asv:an, et.al. is used first 



( 1 ) 10 







(2) from the denominator of a 



o 
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(3) Muller *s method was not used in computing the 

preceding estimate e (It was found that consecutive 
application of this method gave poor convergence.) 

In general, this scheme worked v;ell, ensuring six place 
accuracy within five to eight iterations (depending on the 
accuracy of the initial guess). At first, to ensure uniform con- 
vergence (after the first three arbitrary estimates, ^v^) , 

an estimated ’’eigenvalue" obtained by either method xjas first 
averaged with the most accurate one previously calculated (i.e., 
that one which gave the smallest value of \^^\) before reevaluating 
the stability equations. Hev/ever, since this procedure v;as found 
to increase (by as much as 50%) the num}:>er of iterations required 
for convergence, it was modified to average only xjhen the condition 
of uniform convergence was violated (i.e., when the condition 
|y -V 1 > Iy .1 was not satisfied). Vdien Muller’s method v;as used 
alternately as discussed, it was found to be very^ effective, usually 
refining the "eigenvalue" estimate so as to reduce |y^| by several 
orders of magnitude. Thus, the net result is an iteration scheme 
V7hich ultimately but sporatically converges, often in a nonmontonic 
fashion , 

An example of the convergence is given belov: for = 90®F, 

T - 60^r, M - 0, Q - 50 \ h, = .04, u. = .999, Re = 7000, ex., = 0, 
a_ = .2 2233, o) - 3.h821>:10 ^ (lower branch including density 
f 2 uctuations , r / 0) , 
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Iteration 


"r 


(1) 


|Y 1 

’ o ' 


Guess 1 


.081212293 


.017239455 


.432x10“^ 


Guess 2 


.082024415 


.017411849 


.342x10“^ 


Guess 3 


.081618354 


.017584244 


.139 


1 


.081367993 


.037209441 


.176x10"^ 


2 


.080728449 


.017008234 


.332x10"^ 


3 


.080769324 


.017018352 


.323x10“'^^ 


4 


.080769138 


.017018221 


.311x10“^ 


next estimate. 


.080769140 


.017018221 





table E.S Conversence History of an Eigenvalue Search 



In keeping with previous discussion of obtainable and re 
quircd accuracy, iteration was teriuinatcd when both 




and 



1 



ev 



IR 



ev 



(j+l)R 



and 



\ 






ev 



(j+i)i 



I 



This criterion usually ensured "eigenvalue" accuracy to six places 
and "zero" wall boundary conditions of order (10“^) or less. 
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